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Abstract 

Computational  geometers  invent  algorithms  to  represent  and  manipulate 
geometrical  objects.  The  applications  of  such  algorithms  range  from  the 
creation  of  diagrams  for  documents  to  the  visualization  of  steps  in  the  proofs 
of  correctness  of  motion  planning  algorithms  for  robotics.  An  ideahzation 
of  the  computational  geometric  component  of  each  of  these  applications  is 
the  problem  of  editing  geometric  objects  composed  of  points  and  Hnes  in 
the  plane.  The  editing  algorithms  employed  must  not  lead  to  topologically 
invalid  conclusions,  which  could  result  from  round-off  error  in  munerical 
approximations.  We  model  the  correct  graphical  editing  problem  as  the 
creation,  and  iterative  modification  and  solution,  of  the  set  of  equations 
induced  by  a  network  whose  nodes  are  relations  constraining  the  values  of 
arcs  which  are  labelled  by  portions  of  geometric  objects.  This  survey  focuses 
on  programming  languages  and  systems  which  do  constraint  satisfaction, 
and  on  numerical,  algebraic  and  automated  deduction  methods  which  may 
be  applied  to  solve  the  systems  of  equations  implied  by  constraint  networks. 


This  work  is  partially  supported  by  National  Science  Foundation  grants  DCR-84- 
01898  and  DCR-84-01633,  and  in  part  by  U.S.  Army  Contract  DAAB027-82-K- 
J196. 


Acknowledgement 

I  thank  Professor  Chee  Yap  for  posing  the  problem  of  this  paper  and  direct- 
ing my  reading  in  this  area,  Professors  Bud  Mishra  and  Edmond  Schonberg 
for  critical  reading  of  several  versions  of  this  paper,  and  Professor  Jim  Dem- 
mel  for  informative  conversation  regtirding  numerical  methods.  I  would  also 
like  to  thank  the  Ada  and  Setl  projects  of  Professors  Edmond  Schonberg 
and  Robert  B.K.  Deweir,  for  providing  me  with  a  congenial  working  environ- 
ment and  excellent  criticaJ  opinions  in  the  area  of  programming  language 
design  and  implementation.  Any  remaining  errors  in  this  paper  are  my  own 
fault,  of  course. 


2  CONTENTS 

Contents 

1  Topologically  correct  editing  5 

1.1  The  problem  in  general     5 

1.2  The  geometric  relations  considered 7 

1.2.1  Numbers     8 

1.2.2  Angles 9 

1.2.3  Points 9 

1.2.4  Lines 10 

1.2.5  Frames     12 

1.3  A  canonical  example 12 

1.4  Implementing  a  correct  graphical  editor 16 

2  Constraint  network  manipulations  19 

2.1  Constraint  networks  are  undirected  graphs 19 

2.2  Local  graph  techniques 20 

2.2.1  Propagation  of  known  states     21 

2.2.2  Propagation  of  degrees  of  freedom 24 

2.2.3  Redundant  views 25 

3  Numerical  methods  26 

3.1  Linearizing  constraints 28 

3.1.1  How  to  linearize  constraints 28 

3.1.2  Constraint  linearization  in  SKETCHPAD 28 

3.2  Minimizing  least  squares 30 

3.3  Relaxation 32 

3.3.1  Defining  relaxation 32 

3.3.2  Relaxation  in  SKETCHPAD  and  THINGLAB 33 

3.4  Newton-Raphson  iteration 33 

3.4.1  Details  of  the  method 33 

3.4.2  A  way  of  programming  Newton  iteration 35 

3.4.3  Constraint  systems  which  use  Newton-Raphson  iter- 
ation    36 


CONTENTS  3 

4  Algebraic  methods  37 

4.1  Variations  on  Gaussian  elimination 37 

4.1.1  The  LELAND  Elimination  Algorithm 37 

4.1.2  Algebraic  transformations  in  MaGRITTE     39 

4.1.3  Leler's  augmented  term  rewriting  system 40 

4.2  Advanced  algebraic  methods     41 

4.2.1  van  Wyk's  observations  on  algebraic  complexity     .  .  41 

4.2.2  Grobner  bases  and  ideal-theoretic  methods 42 

4.2.3  Algebraic  nimabers 46 

5  Finite-precision  arithmetic  and  graphics  49 

6  Constraint  languages  50 

6.1  The  SKETCHPAD  graphical  editor 50 

6.2  A  proposal  of  Wilkes 52 

6.3  Constraint  languages  and  applications  at  MIT 53 

6.4  The  THINGLAB  graphical  simulation  language 54 

6.5  The  MAGRITTE  graphical  editor 56 

6.6  The  Ideal  typesetting  language 56 

6.7  The  Juno  graphical  editor     57 

6.7.1  Elements  of  the  language 57 

6.7.2  Juno  code  for  an  equilateral  triangle     58 

6.8  The  BERTRAND  production  system  language 58 

6.8.1  Rules  in  BERTRAND 58 

6.8.2  Types  in  BERTRAND 59 

6.8.3  Evaluation  order  in  BERTRAND 60 

7  Conclusion  61 
References  62 


LIST  OF  TABLES 


List  of  Figures 


1  Topological  error  caused  by  approximation 6 

2  Line  types 11 

3  The  intersection  of  two  lines 13 

4  Addition  and  subtraction 14 

5  A  +  B  =  C*  D 14 

6  The  intersection  of  two  lines,  in  detail 15 

7  The  constraint  x  +  x  —  4 23 

8  Local  propagation  fails  here 24 

9  Local  propagation  succeeds  here     26 

List  of  Tables 

1        Graphical  editors  and  languages  and  their  methods 17 


1      Topologically  correct  editing 
1.1      The  problem  in  general 

Current  computer  graphics  systems  use  numerical  approximations  when 
calculating  the  position  of  objects  in  relation  to  each  other.  These  ap- 
proximations may  result  in  topologically  invalid  conclusions,  when  the  \in- 
certainty  induced  by  round-off  error  is  not  accounted  for.  Suppose,  for 
example,  that  we  define  three  lines 

Lx-y   =    X 
Lb'-V    =    -x  +  2 
Lc-.y    =    1.1 

and  define  a  certain  point  px  as  the  intersection  of  La  and  Lb-  By  algebraic 
solution  of  the  linear  equation  system  we  see  that  p^  lies  below  (1,1.1).  But 
iising  numerical  methods  with  fixed-precision  arithmetic,  we  niay  instead 
conclude  that  px  is  above.  Since  the  algebraic  solution  is  niore  accurate,  we 
see  that  simple  numerical  approximation  results  in  an  error  in  computing 
the  discrete  relationship  {below,  above}  between  two  objects.  The  situation 
is  illustrated  in  Figure  1. 

The  problem  is  to  create  a  computer  program  to  edit  point  and  line 
objects  and  groups  of  objects  in  the  plane,  in  a  correct  fashion  (to  be  defined 
in  a  moment). 

Objects  are  partially  specified  by  a  set  of  geometric  relations  between 
objects.  This  set  may  be  viewed  as  a  network  of  constraints  on  the  objects. 
In  a  constraint  network,  the  arcs  are  labels  of  portions  of  objects,  and  the 
nodes  are  geometric  relations. 

An  editing  step  is  performed  by  modifying  the  constraint  network,  then 
computing  a  satisfying  assignment  of  values  to  variables  in  the  system  of 
equations  induced  by  the  network.  Depending  on  how  much  information  is 
"filled  in",  the  result  of  an  editing  step  is  either  a  reduced  set  of  equations, 
with  some  unknowns  remaining,  or  a  complete  set  of  satisfying  assignments 
to  all  va.riables  occurring  in  the  equations.  As  a  refinement,  we  recursively 
apply  the  process  to  isolatable  subnetworks. 
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a.  The  actual  intersection  point  is  at  (1,1) 


b.  Shift  due  to  approximation  of  intersection  point 


Figure  1:  Topological  error  caused  by  approximation. 
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By  correct  we  mean  that  no  network-satisfying  assignment  leads  us  to 
conclude  that  a  discrete  relationship  holds  between  two  objects,  when  we 
would  conclude  that  the  relationship  did  not  hold  if  we  performed  an  exact 
algebraic  solution  of  the  system. 

The  problem  above  was  suggested  to  me  by  Prof.  Chee-Keng  Yap,  in  the 
course  of  an  independent  study  on  computer  algebra  (Spring,  1986).  During 
this  time  we  worked  on  a  language  for  describing  constrained  geometric 
objects. 

1.2      The  geometric  relations  considered 

The  complexity  of  correct  editing  increases  with  the  number  of  dimensions^, 
and  is  unsolvable  for  unrestricted  arithmetic  relationships  on  the  domain 
of  integers  [8].  We  restrict  ourselves  to  the  plane  and,  most  importantly, 
we  restrict  our  investigation  to  methods  which  apply  to  a  particular  set  of 
geometric  relationships,  which  we  outline  below. 

The  following  objects  and  relational  operators  are  those  expressed  in 
Ericson  and  Yap's  propK)sed  LINEToOL  language  [ll]. 

Picture  a  plane  on  which  a  collection  of  Unes  and  points  is  displayed. 
These  lines  and  points  are  organized  into  movable  figures,  which  we  call 
frames.  Within  a  frame,  lines  and  points  are  confi^red  in  relation  to  each 
other  by  geometric  assertions,  such  as  the  statement  that  point  p  is  formed 
by  the  intersection  of  lines  Li  and  Lj.  Frames  are  projected  onto  the  plane, 
where  the  projection  is  defined  by  a  translation  point  and  a  rotation  angle. 
The  translation  point  locates  the  origin  of  the  virtual  plane  defined  by 
the  frame  on  the  absolute  plane.  The  rotation  angle  gives  the  angle  of 
the  frame  plane  i-aocis  with  respect  to  the  absolute  plane  z-axis  and  the 
absolute  origin. 

We  may  think  of  many  different  ways  of  describing  a  certain  type  of 
object.  For  example,  a  line  may  be  defined  by  a  linear  equation  y  =  mx  +  b, 
or  it  may  be  defined  by  two  points  through  which  it  passes,  or  it  may  be 
defined  by  a  point  and  an  angle;  there  are  yet  other  descriptions.  But  for 


^1  don't  want  to  say  much  more  than  that,  because  it  is  basically  a  qualitative  statement 
which  I  can't  give  general  support  for,  but  is  used  to  justify  restricting  the  problem  to 
the  plane. 
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purposes  of  definition,  we  choose  a  canonical  representation  for  each  type 
of  object  we  will  consider.  Alternate  representations  may  then  be  defined 
by  equations  involving  the  canonical  representaitons. 

1.2.1     Numbers 

The  most  primitive  type  of  object  required  in  the  specification  of  geometric 
objects  is  the  number.  There  are  various  subtypes  of  number  together  with 
their  representations,  such  as  real,  complex,  rational  and  algebraic  numbers. 
For  the  purpose  of  specifying  objects  by  their  relation  to  each  other,  it  does 
not  concern  us  how  a  number  is  represented:  for  topological  correctness 
to  hold,  the  statement  of  relationships  between  geometric  objects  such  as 
points  and  lines  should  be  independent  of  the  representation  chosen  for  a 
nimiber  at  a  particular  point  in  the  computation.  We  also  generally  do  not 
care  about  which  subtype  of  number  is  required;  the  "computing  engine" 
should  choose  the  correct  subtype  in  the  course  of  computing  relationships. 
The  following  are  some  number  representations  that  may  be  useful. 

•  Fixed- precision  floating-point  numbers,  represented  by  a  word  or  two 
of  computer  memory,  together  with  non-cissociative  rounding  or  trun- 
cating machine  arithmetic. 

•  Rational  numbers,  represented  by  pairs  of  integers  which  form  a  quo- 
tient, together  with  greatest  common  divisor  nonnalizing  rational 
arithmetic. 

•  Algebraic  numbers,  represented  by  a  defining  polynomial  with  rational 
or  integer  coefficients  plus  an  isolating  interval  for  a  particular  root 
of  the  polynomial  represented  by  a  pair  of  rationals.  The  topic  of 
algebraic  number  arithmetic  is  quite  complex  and  will  be  elaborated 
in  a  section  below. 

We  assimae  the  presence  of  the  relations  —,  ^,  and  <  on  numbers, 
as  well  SiS  the  relations  of  addition,  multiplication  and  exponentiation.  For 
example,  we  view  addition  on  (a,  6,  c)  as  the  relation  a  +  6  =  c.  The  idea  here 
is  that  we  are  specifying  geometric  objects  indirectly  by  their  relationship 
to  each  other.    The  arithmetic  operations  plus,  times  and  exponent  are 
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viewed  as  relations  rather  than  functions  because  it  is  more  convenient.  A 
function  is  a  one-directional  operation  which  holds  when  we  have  values  for 
all  of  its  arguments.  Thus  Add(i,y)  =  i  +  y  is  a  function  which  takes  on 
a  value  when  we  know  x  and  y.  The  relation  of  addition  x+  y  =  z  is  much 
more  flexible:  it  tells  us  all  when  we  know  i  and  y,  or  i  and  z,  or  y  and 
2,  and  it  tells  us  something  even  when  we  know  only  one  of  z,  y  or  z,  or  if 
we  have  only  partial  information  about  some  subexpression  of  x,  y  or  z,  in 
the  case  that  these  are  expressions  rather  than  simple  variables. 

1.2.2      Angles 

ATigles  'are  simply  numbers,  together  with  arithmetic  modulo  the  period 
of  the  chosen  units  {2n  in  the  case  of  radians,  360  in  the  case  of  degrees). 
Again,  for  the  purpose  of  surveying  the  literature,  we  do  not  care  about 
the  exact  representation. 

Between  numbers  and  angles  we  define  a  coercion  relation 

RNxi^iO,)  =  »■  =  number(a)  =  a  =  angle(r) 

where  r  is  a  real  number  and  a  is  an  angle. 

Between  numbers  and  angles  we  have  the  cosine  relation 

Rcot{r,<i)  =  r  =  cos  a  =  a  —  arccosr 

A  special  angle  is  tt: 


Rf{a)  =  a  =  ir 


1.2.3     Points 


A  point  I  is  defined  by  two  nimabers  (11,12). 

The  relation  of  Euclidean  distance  d  between  two  points  x  =  (xi.Zz) 
and  y  =  (yi,y2)  is  defined  by 

REucUd{x,y,d)  =  d=  [(ii  -  I2)*  +  (yi  -  y:)')]^ 

We  assimae  the  relation  of  distance  between  a  point  and  a  line. 
We  also  assume  the  relation  of  intersection,  that  is,  of  defining  a  point 
as  the  intersection  of  two  lines.  To  define  this,  it  is  more  convenient  to  use 
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the  homogentous  eoordinatt  techniques  described  by  Newman  and  Sproull 
[28].'  Let  MS  then  represent  a  point  p  =  (x,  y)  by  the  row  vector  v  = 
[wx  wy  w],  where  w  is  called  a  scale  factor.*  Let  us  represent  the  line 
with  equation  ax  +  by  +  c  =  0  by  the  column  vector 


1  = 


a 
b 
c 


Then  the  relation  defining  that  a  point  v  to  be  the  intersection  of  two  lines 
7,  A  is  simply  (where  x  is  the  vector  cross  product): 

■Rx(v,'r,A)  =  V  =  7  X  A 

1.2.4     Lines 

A  line  L  is  a  4-tuple  {lp,ltyp€,rp,rtype)  where  Itype  and  rtype  are  the  type 
of  the  left  and  right  points  Ip  and  rp,  and  a  type  is  one  of 

{open,  closed,  through} 

Open  and  closed  signify  endpoints.  A  closed  endpoint  is  one  that  the  line 
terminates  on.  An  open  endpoint  is  one  that  the  line  terminates  just  before. 
A  through  point  is  a  point  that  the  line  passes  through.  Figure  2  displays 
the  six  different  kinds  of  line  object  we  may  represent  in  this  fashion. 

On  lines  we  assume  relations  of  slope,  intersection  (see  above),  distance 
to  a  point,  left  and  right  endpoint  location  and  reversal.  In  case  we  manage 
to  specify  that  two  parallel  lines  intersect,  in  the  process  of  relating  a  set 
of  geometric  objects,  then  two  courses  of  action  are  open.  If  there  is  more 
than  one  possible  "completion"  of  the  scene  given  the  initial  specification, 
then  choose  a  completion  which  does  not  cause  parallel  lines  to  intersect. 
Otherwise  we  conclude  that  the  description  is  inconsistent. 

Since  there  are  many  way?  of  specifying,  say,  a  line,  it  is  not  clear  what  the  best  internal 
representation  of  geometric  objects  might  be.  However,  matrix  equations  on  homoge- 
nous coordinates,  which  are  standard  in  computer  gfraphics,  may  be  the  best  internal 
representation  of  geometric  object  specifications  within  a  graphical  editor. 

One  practical  significance  of  the  scale  factor  is  that  it  allows  us  to  normalize  steps  fixed- 
precision  computations  to  avoid  round-off  or  out-of-scale  errors.  We  do  not  concern 
ourselves  with  this  here. 


1.2     The  geometric  relations  considered 


11 


(4,  4),  through 


(1,  1),  through 
a.  Full  line 
(4,  4),  closed 


(4,  4),  open 


(1,  1),  through  y*  (1.  l).  through 


b.  Half  lines 


(4,  4),  closed  «  (4,  4),  open  «  (4,  4),  open 


[1,  1),  closed 


''(1,  1),  closed  c/(l,  1),  open 

c.  Line  segments 


Figure  2:  Line  types. 
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1.2.5     Frames 

A  frame  F  =  {S,  p,  a),  where  5  is  a  set  of  lines  and  points,  and  the  projection 
of  the  frame  onto  the  plane  is  defined  by  the  point  p  and  angle  a,  which  set 
the  translation  from  the  origin  and  rotation  of  coordinates  of  objects  in  S. 
The  translation  point  (i,y)  fixes  the  location  of  the  frame  origin  (0,0)  with 
respect  to  the  plane  origin  (0,0).  The  rotation  angle  0  fixes  the  angle  of 
the  X  —  y  axis  of  the  frame  plane  with  respect  to  the  absolute  plane,  where 
the  absolute  plane's  x-axis  lies  along  the  0-degree  radial  and  the  absolute 
plane's  y-axis  lies  along  the  90-degree  radial,  and  the  frame  plane's  i-axis 
lies  along  the  fl-degree  raxlial,  and  the  frame  plane's  y-axis  lies  along  the 
6  -\-  90-degree  radial.  Frames  may  be  related  to  each  other  by  assertions 
relating  their  translations  and  rotations. 

On  frames  we  assume  relations  of  translation,  transformation  and  merg- 
ing. 

1.3      A  canonical  example 

We  will  use  the  problem  of  defining  a  point  p  =  (x,  y)  as  the  intersection  of 
two  lines  Ly^  and  Lb  as  a  running  example  in  this  paper.  Lyt  and  Lb  are 
defijied  by  the  linear  system 

OiX  -I-  Ojy  +  as    =    0 
bix  +  bty+bs    =    0 

We  think  of  the  intersection  point  (x,y)  as  being  constrained  by  the 
defining  equations  of  the  two  lines.  Similarly,  if  the  intersection  point  is 
fixed,  but  some  coefficients  of  the  line  equations  are  not,  then  we  thiak  of 
the  known  intersection  point  as  constraining  the  line  equations. 

Figure  3  says  that  point  Px  is  defined  by  the  intersection  of  two  lines  L^ 
and  Lb,  where  intersection  is  defined  by  the  relation  R^{Lx,Lb,Px)  and 
Rx  is  defined  as  follows: 


Rx{aix^a2y+ai,bix+b2y+bs,{xi,yi))  =  Xi  =  — — Ayi  =  — 


1.3     A  canonic ai  exampJe 
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Figure  3:  The  intersection  of  two  lines. 

Note  that  it  is  often  possible  to  trivially  express  a  number  of  relations 
with  one  relation.  For  example,  Figure  4  illustrates  how  to  apply  the  addi- 
tion relation 

R+{xi,X2,Xi)  =  Xi  +  X2  =  Xi 

to  (a,  6,  c)  to  express  a  +  b  =  c  and  to  (6,  c,  a)  to  express  a  —  b  =  c. 

We  do  not  need  to  define  i2=.  Equivalence  in  a  constraint  network  can 
be  elided  by  merging  output  arcs.  For  example,  the  equation  a  -\-  b  =  c  *  d 
is  shown  in  Figure  5. 

Figure  6  shows  the  fully-instantiated  constraint  network  corresponding 
to  Figure  3.  This  network  is  derived  from  the  equations  of  the  solution  for 
I  and  y. 

An  instance  of  the  intersection  problem,  if  we  suppose  that  i  =  12, 
ai  =  3,  Cs  =  4  and  fej  =  1,  is  to  take  the  constraint  network  of  Figure  6 
and  apply  the  relation  Rx  to  obtain  the  following  system  of  equations: 


12 

y 


0263  —  1462 
02^1  —  3 

146i  -  363 
02^1  ~  3 


and  find  simplified  equations  for  the  remaining  unknowns  aj,  61  and  63.  For 
example,  for  Cj  we  have 

-        22 
°'~  1261  -  63 
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A 
B 


B 
C 


a.  A  +  B  =  C 


b.  A-  B  =  C 

Figure  4:  Addition  and  subtraction. 


A 
B 


C 
D 


Figure  5:  A  +  B  =  C  *  D. 


1.3     A  c&DonicaJ  examp/e 
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oi  aj  as  61  bih  X  V 


+ 

1 

+ 

+ 

+ 

Figure  6:  The  intersection  of  two  lines,  in  detail. 
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1.4     Implementing  a  correct  graphical  editor 

We  break  down  the  issues  of  implementing  a  constraint-based  graphical 
editor  as  follows.  The  paper  will  discuss: 

•  Methods  of  reducing  and  transforming  the  constraint  graph  to  satisfy 
its  component  relations. 

•  Numerical  methods  for  approximately  solving  constraint  networks. 

•  Algebraic  methods  for  operating  on  the  finite  system  of  equations 
induced  by  a  constraint  network,  to  reduce  the  system  to  a  simpler 
one,  to  find  zeroes  of  its  constituent  polynomials,  and  to  compute 
relationships  between  geometric  objects.  Included  under  this  heading 
are  elimination  theory,  computation  of  Grobner  bases  for  ideals  of 
systems,  term  rewriting  systems,  and  algebraic  number  arithmetic. 

•  Considerations  from  computer  graphics  for  the  correct  or  at  least  non- 
deceptive  display  of  high-accuracy  graphical  information  on  finite- 
accuracy  displays. 

•  Design  of  the  graphical  object  description  language,  and  association  of 
constraint-network  semantics  to  elements  of  the  description  language. 

Table  1  summarizes  the  methods  in  use  by  various  graphical  editors  and 
languages. 

The  network  problem  is  obviously  an  isomorph  of  the  equational  prob- 
lem. The  "network  manipulation  approach"  is  a  school  of  thought  which 
hais  been  developed  in  the  following  works:  Sutherland  [36];  Borning  [3]; 
Gosling  [13];  and  Sussman,  Steele  and  Zippel  [33,35,44]. 

Numerical  methods  have  the  failing  that  that  they  are,  in  the  usual 
implementation  in  which  only  an  absolute  level  of  precision  is  prescribed 
for  a  calculation,  not  linear-order  preserving.  By  linear  order  preserving  we 
mean  that  no  sequences  of  expression  evaluations  causes  us  to  conclude  that 
a  <  b  when  in  fact  a  >  b.  These  methods  are  fast  (but  slower  than  network- 
manipulation  methods,  when  they  work).  Lazard  notes  that  numerical 
methods  do  not  give  the  number  of  roots  of  a  system  of  equations,  nor  do 
they  work  on  extensions  of  finite  fields  [19].  An  example  of  the  fallibility  of 
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System 

Year 

Author 

Methods 

Sketchpad 

Sutherland 

1963 

Least  mean  squares  fit  relaxation 
Network  manipulation 

Thinglab 

Borning 

1979 

Relaxation 

Network  manipulation 

Ideal 

van  Wyk 

1980 

Gaussian  elimination 

Constraints 

Steele 

1980 

Network  manipulation 

Magritte 

Gosling 

1983 

Gaussian  elimination 

Juno 

Nelson 

1985 

Newton-Raphson  iteration 

Table  1:  Graphical  editors  and  languages  and  their  methods 


numerical  methods  is  as  follows.  The  JUNO  graphical  editor  [27],  which  uses 
the  Newton-Raphson  method  to  solve  constraints,  allows  one  to  constrain 
two  line  segments  to  be  parallel.  It  also  allows  one  to  declare  the  length 
or  relative  length  of  a  line  segment.  Nelson  notes  that  "if  the  user  makes 
a  long  segment  parallel  to  a  short  segment,  the  solver  tends  to  collapse  the 
short  segment  to  a  point."  In  other  words,  the  length  of  the  shorter  line 
segment  is  set  to  zero. 

Language  design  is  largely  a  matter  of  taste,  because  the  limiting  factor 
in  the  solution  of  our  main  problem  as  posed  in  the  first  section  is  the 
availability  of  fast  algorithms  for  solving  systems  of  equations.  However, 
we  will  discuss  several  distinct  and  innovative  language  designs. 

The  equational  approach  has  a  much  richer  theoretical  background,  but 
only  recently  has  there  been  much  pragmatic  success  in  topologically  correct 
methods. 

Elimination  methods  eliminate  successive  variables  un  til  one  equation 
in  one  variable  remains.  All  roots  are  found.  Such  methods  can  be  slow 
(though  specific  methods  for  linear  systems  are  fast).  As  Lazard  notes,  "if 
the  system  consists  of  equations  of  degree  rf  in  n  variables,  the  final  equa- 
tion has,  in  general,  degree  cP"  ."  In  other  words,  in  general  the  obvious 
method  of  ehmination  is  intractable.  We  will  group  all  such  methods  under 
the  heading  "Gaussian  elimination",  though  the  term  has  a  more  restrictive 
technical  definition. 
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The  family  of  Grobner  basis  methods  provides  a  conceptually  differ- 
ent approach  to  elimination  which  is  much  faster  under  certain  conditions. 
These  methods  can  be  shown  to  be  similar  to  elimination  [20],  but  the 
problem  is  phrased  differently.  The  complexity  of  Grobner  basis  methods 
in  general  is  an  open  question  [26]. 

Term  rewriting  systems  are  simply  a  family  of  methods  for  applying 
equational  axioms  to  initial  sets  of  hypotheses.  This  paradigm  is  employed 
by  Leler  to  specify  constraint  satisfaction  systems  [21]. 

An  algebraic  number  is  defined  a  particular  zero  of  a  given  polynomial.  It 
is  possible  to  do  linear  order  preserving  arithmetic  on  algebraic  numbers. 
With  algebraic  numbers  this  is  done  by  increasing  the  precision  of  the 
intervals  when  necessary  in  the  course  of  performing  arithmetic  operations. 
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2      Constraint  network  manipulations 

In  this  section  we  discuss  a  few  non-numerical  constraint  network  tech- 
niques that  have  been  developed  by  those  authors  whose  systems  use  the 
constraint  graph  as  the  preferred  object  of  computation. 

We  begin  with  a  formal  definition  of  the  ba^ic  constraint  graph  definition 
assumed  by  all  authors.  Then  we  give  names  and  sketches  of  the  algorithms 
for  the  network  techniques.  The  purpose  of  the  techniques  is  to  find  a 
satisfying  assignment  for  the  variables  constrained  by  the  network. 

Successive  summaries  exist  of  the  techniques  we  shall  describe.  The 
most  recent  summary  is  that  of  William  Leler  [21,  chapters  1,2  and  5]. 
Prior  authors  who  have  had  to  summarize  these  techniques  as  part  of  their 
preliminary  thesis  work  include  Borning  [3],  Steele  [33]  and  Gosling  [13].  ■* 


2.1      Constraint  networks  are  undirected  graphs 

Define  a  constraint  relation  R^  to  be  a  set  of  n-tuples  x,  where  n  is  the 
arity  of  the  relation,  and  x  G  JS".  (Of  course,  the  set  may  be  infinite.)  Each 
constraint  relation  R^  has  associated  with  it  a  defining  algebraic  expression 
^(x),  such  that 

<^(x)  =  0  =  X  G  /2^ 

We  leave  unspecified  the  exact  class  of  expressions  of  which  ^(x)  is  a 
member,  except  to  say  that  we  intend  the  class  that  includes  multivariate 
polynomials  and  involves  trigonometric  and  logarithmic  functions  (but  no 
differentials  or  integrals). 

Note  also  that  R^,  and  Rz-\  are  relations  of  arity  1  denoting  the  con- 
stants 0  and  1.  That  is,  by  the  above  definition,  z  G  72,  =  a;  =  0,  and 
I  6  Rx-\  =  X  —  1  =  0,  that  is,  x  6  Rx-i  =  i  =  1.  Also  R^  where 
(/)(x)  =  Xi  +  Zj  —  Xs  denotes  the  addition  relation. 

Define  a  constraint  graph  to  be  an  undirected  graph  G  =  (V,  £■,  A,cr,  L) 
where 


*Since  Sutherland  wa«  the  first  to  apply  such  techniques,  he  did  not  have  to  Bummarize 
anything. 
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1.  The  vertices  V  denote  constraint  relations,  as  induced  by  the  labelling 

X:V  —*{<i>„4>2,..-4>n} 
where  the  <f>i  are  defining  algebraic  expressions. 

2.  The  edges  E  denote  variables,  such  that  the  partial  assignment  map 

a:  E  — y  R. 
gives  the  current  aissigned  value  of  the  variable  denoted  by  the  edge. 

3.  For  every  vertex  t;  G  V  its  set  of  incident  edges  e  is  ordered  Ci,  cj, ...,  e^ 
where  k  is  the  arity  of  A(v).  Also,  let  Li{v)  denote  Ci  in  this  ordering. 

4.  It  may  be  possible  for  edges  E  to  have  one  or  both  of  their  endpoints 
dis connect td.  An  edge  with  one  endpoint  disconnected  is  called  a 
free  variable.  An  edge  with  both  endpoints  disconnected  is  called  a 
redundant  variable. 

5.  A  completion  of  a  partial  assignment  map  cr  is  a  total  assignment 
map  o'  that  extends  o.  A  partial  assignment  map  a  is  consistent 
when  there  exists  a  completion  a'  of  o  such  that  for  all  v  G  V  the 
following  holds: 

Rxi.){o'{Li{v)),a'{L,{v)),...,a'{Lr,{v))) 

2.2      Local  graph  techniques 

The  following  techniques  are  used  to  compute  portions  of  satisfying  cissign- 
ments,  which  have  the  property  that  they  operate  on  information  in  the 
neighborhood  of  individual  vertices.  Note  that  the  neighborhood  of  a  given 
vertex  is  just  that  set  of  vertices  with  edges  to  the  given  vertex. 

•  Propagation  of  known  states,  from  Boming,  also  called  one-shot  prop- 
agation by  Gosling  at  one  point  and  unplanned  firing  at  another. 

•  Propagation  of  degrees  of  freedom,  so-called  by  Borning,  first  used  by 
Sutherland,  who  called  it  the  one-pass  method  (we  adopt  the  former 
term). 
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•  Mtiltiple  views,  from  Boming,  called  alternate  views  by  Gosling  and 
redundant  views  by  Leler  (we  adopt  the  latter  term). 

2.2.1     Propagation  of  known  states 

Propagation  of  known  states  applies  when  the  values  of  a  sufficient  number 
of  edges  of  a  relational  vertex  have  been  assigned,  so  that  the  values  of 
additional  arcs  may  be  fixed  by  the  definition  of  the  relation.  For  example, 
if  a  vertex  denotes  the  addition  a  +  6  =  c,  and  the  edges  corresponding  to 
variables  a  and  c  are  known  to  be  2  and  4,  respectively,  then  the  edge  for  b 
may  be  assigned  the  value  2,  and  the  vertex  may  be  removed  from  further 
consideration.  In  a  graph  with  enough  pre-assigned  values  for  edges,  with 
vertex  relations  which  are  sufficiently  solvable,  iterating  this  process  over 
the  remaining  vertices  in  the  graph  will  be  sufficient  to  reduce  the  graph 
to  a  single  vertex  whose  partial  assignment  is  satisfying.  This  leads  to  the 
following  algorithm  (compare  Gosling's  presentation  [13,  page  23]). 

Algorithm  LocaiPropagationSolve  {G  =  {V,E,X,<7,L)) 

Input  A  constraint  graph  G. 

Output  G'  =  (F,£^,A,a,  L),  where  a*  is  a  maximal  extension  of  a,  or 
failiu-e  if  a  is  inconsistent  or  G  is  unsatisfiable.  There  may  be  more 
than  one  maximal  extension.  If  so,  the  algorithm  will  return  the 
first  one  found,  but  the  particular  one  foimd  may  be  a  function  of 
extraneous  factors  such  as  the  order  of  storage  of  vertices  in  the  graph 
representation. 

Local  variables  V,  a  set  of  constraint  relation  vertices.  E',  a  set  of  edges. 
a',  a  i>artial  assignment,  v,  a  particular  constraint  vertex  under  ex- 
amination. 

Method 

[Initialise]  a'  :=  a;  E'  :=  E;  V  :=  V 

[E  empty?]  if  i;  =  0  then  return  a" 

[V  empty?]  If  V  =  0  then  return  a"  noting  that  e  E  E  are  redundant 
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[Inconsistent?]  if  3t;  6  V  such  that  v  is  inconsistent  under  </  then 

fail 
[Reduce  E]  if  a"  forces  any  e  e  E  under  the  relations  A  of  F  then 

1.  a*  :=  ty'U  newly  forced  values 

2.  £'  :=  E'—  forced  edges 

3.  goto  [E  empty?] 

[Reduce  V]  if  3v  6  V  such  that  all  edges  of  v  have  values  in  a'  then 

1.  v'~V'-{v} 

2.  goto  [E  empty?] 

Two  constraint  languages,  the  family  of  languages  investigated  by  Steele 
[33],  and  the  THINGLAB  system  of  Borning  [3],  propose  that  one  enumer- 
ate, for  each  possible  case  in  which  the  partiaJ  assignment  determines  the 
relation  uniquely,  a  rule  for  computing  the  remaining  components.  For 
example,  the  following  is  Steele's  language  for  the  multiplication  relation:^ 

(define-primitive-constraint  multiplier  (a  b  c) 
((a)        (and  (is-zero  a)  (set  c  0))) 
((b)        (and  (is-zefo  b)  (set  c  0))) 
((a  b)     (set  c  (times  a  b))) 
((a  c)     (and    (not  (is-zero  a)) 

(is-zero  (remainder  c  a)) 

(set  b  (quotient  c  a)))) 
((b  c)     (and    (not  (is-zero  b)) 

(is-zero  (remainder  c  b)) 

(set  a  (quotient  c  b))))) 

What  this  says  is  that  we  are  defining  a  type  of  constraint  relation  called 
a  multiplier,  with  a  parameter  with  three  components  a,  b  and  c.  If  a  is  in 
the  partial  assignment  of  a  relation  instance  of  multiplier,  and  o  =  0,  then 
add  (c,0)  to  the  partial  assignment.  Likewise  if  6  =  0.  If  a  and  b  are  in  the 
partial  tissignment,  add  (c,a6)  to  the  partial  assignment.  If  a  and  c  are  in 
the  partial  assignment,  a  /  0  and  c  divides  a,  then  add  (6,  ^)  to  the  partial 
assignment.  Likewise  if  6  and  c  are  in  the  partial  assignment. 


^I  have  expanded  the  function  names  for  thoee  not  familiar  with  Lisp  Machine  LiSP,  but 
otherwise  the  code  is  as  presented  in  Steele's  thesis  [33]. 
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Figure  7:  The  constraint  x  +  x  =  4. 

We  note  pragnaatically  that  if  there  are  N  components  to  the  parameter 
X  of  the  relation  i?  of  a  vertex,  then  there  will  be 


=  2^^  -  1  (1) 


possible  combinations  of  specified  components  in  non-total  partial  assign- 
ments. Hence  there  are  that  many  different  ways  that  a  relation  may  be 
underdetermined,  in  general.  Clearly  a  language,  such  as  Steele's  shown 
above,  which  proposes  that  the  iiser  make  this  enumeration,  will  lead  to 
very  large  programs. 

Both  because  of  this  difficulty,  and  because  Steele  and  Boming's  lan- 
guages do  not  have  a  way  of  testing  whether  two  edges  actually  denote  the 
same  unknown  (which  is  possible  if  two  edges  have  the  same  source),  Steele 
and  Borning  restrict  themselves,  in  their  application  of  local  propagation, 
to  reducing  nodes  for  which  the  number  of  knowTi  values  is  n  —  1  where 
n  is  the  arity  of  the  node.  For  example,  Leler  gives  the  example  of  the 
constraint  that  x  +  x  =  4,  depicted  in  Figure  7.  Clearly  this  is  solvable, 
though  the  number  of  known  values  is  n  —  2,  if  we  note  the  identity  of  the 
unknowns. 

For  situations  best  described  by  a  non-collapsing  system  of  simultaneous 
linear  equations,  the  next  most  powerful  technique  prescribed  by  Borning 
and  discussed  by  Gosling  and  Leler  is  that  of  propagation  of  degrees  of  free- 
dom, discussed  in  the  next  section.  In  particular,  Gosling  gives  an  example 
of  a  constraint  network  for  which  local  propagation  (with  or  without  testing 
for  the  equality  of  unknowns,  and  with  or  without  solutions  for  numbers  of 
known  values  less  than  n  —  l)  does  not  work.  The  example  corresponds  to 
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+ 

B 

+ 

10 


Figure  8:  Local  propagation  fails  here. 

the  simultaneous  linear  eqiiation  system 

5  +  r    =  B 
B  +  T    =10 

This  example  is  depicted  in  Figure  8. 


2.2.2     Propagation  of  degrees  of  freedom 

Recall  that  a  free  variable  in  a  constraint  graph  is  an  edge  with  one  endpoint 
connected  to  no  vertex.  The  number  of  degrees  of  freedom  in  a  constraint 
graph  is  simply  the  number  of  free  variables.  The  idea  of  the  method 
of  propagation  of  degrees  of  freedom,  as  variously  described  by  Boming, 
Gosling  and  Leler,  appears  to  be  to 

1.  Remove  a  vertex  v  from  the  graph  G  which  is  solvable  modulo  a  free 
variable,  that  is,  for  which  forced  assignments  to  its  edges  would  be 
computable  if  a  free  variable  incident  on  the  vertex  were  aissigned  a 
value,  and  let  Gs  be  the  resulting  graph  If  no  such  vertex  exists,  then 
fail. 

2.  Attempt  to  solve  the  resulting  subgraph  Gs,  first  by  the  method  of 
propagation  of  known  states,  then  recursively  by  this  method,  then 
by  relajtation  or  Newton-Raphson  iteration,  otherwise  fail. 

3.  Use  the  resulting  solution  for  the  subgraph  Gs  to  solve  v  and  hence 
solve  G. 
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Sketchpad  mbcs  the  above  method,  where  step  2  invokes  relaxation  if 
the  above  method  fadls  recursively,  and  step  3  invokes  propagation  of  known 
states  [13,  pages  38-39]. 

Upon  reading  Gosling,  Leler  and  Boming's  accounts,  the  author  was  un- 
able to  concoct  an  example  applying  the  above  method,  which  would  either 
succeed  where  propagation  of  known  states  failed,  or  succeed  where  relax- 
ation, Newton  iteration  or  Gaussian  elimination  wa*  otherwise  required. 
The  cited  accounts  either  assumed  that  propagation  of  known  states  would 
fail  on  the  example  of  Figure  7,  or  that  propagation  of  degrees  of  free- 
dom was  being  used  as  a  compilation  step,  for  example  to  compile  the 
expression  C  =  {F  —  32)/ 1.8)  from  the  constraint  graph  for  the  expression 
1.8C  +  32  =  F. 

2.2.3     Redundant  views 

If  the  constraint  solver  is  not  powerful  enough  to  satisfy  a  given  set  of 
constraints  as  specified,  it  may  be  necessary  for  the  programmer  to  over- 
constrain  (from  a  mathematical  point  of  view)  the  problem  to  the  extent 
necessary  to  revive  the  solver.  This  is  classified  as  a  local  technique  in  the 
sense  that  it  is  a  meta-technique  (also  called  a  hack)  applied  by  the  user  of 
a  constraint-satisfaction  system  such  as  THINGLAB  [3]  or  Steele's  languages 
[33]  in  order  to  revive  that  portion  of  the  constraint  solver  which  does  local 
propagation  of  known  states. 

To  extend  Gosling's  example  of  Figure  8,  what  we  wo\ild  do  is  overcon- 

strain  the  system 

5  +  T    =  B 

B  +  T    =10 


by  adding  the  equations 


10    =5  +  T2 
T2    =  2  *  r 


See  Figure  9  for  the  resulting  network.  This  addition  causes  the  method  of 
local  propagation  of  known  states  to  succeed. 


26 


3    NUMERICAL  METHODS 


Figure  9:  Local  propagation  succeeds  here 


3      Numerical  methods 

In  general,  the  problem  of  constraint  satisfaction  amounts  to  translating  a 
constraint  network  into  a  system  of  M  equations  in  N  unknowns,  and  then 
doing  something  with  the  resulting  system.  If 

•  M  >  N,  then  the  system  is  overconstrained,  and  we  may  choose  to 
either  declare  the  system  inconsistent,  or  find  an  assignment  to  the 
variables  which  represents  an  interpolation  of  the  constraints. 

•  M  =  N,  then  we  seek  an  approximation  to  the  zeroes  which  is  exact 
in  some  sense. 

•  M  <  N,  then  perhaps  we  can  reduce  M  even  further.  But  at  any  rate, 
what  we  have  is  a  system  of  constraints  on  an  infinite  set  of  solutions, 
with  no  distinct  solution.  We  may  choose  for  application  purposes 
to  pick  a  "pleasing"  solution  from  this  set  using  some  modification  of 
the  techniques  discussed  below,  but  otherwise  nothing  further  should 
be  done. 

So  from  the  nimierical  viewpoint,  the  problem  is  to  solve  or  optimize 
a  system  of  equations,  that  is,  to  find  assignments  to  the  Xj  such  that  for 
1  <  »  <  M, 

/,{zi,X2...iAr)  =  0, 
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or  more  compaw;tly,  to  solve  for  x  the  equation 

f  (x)  =  0 

where 

X     =      [z,       Xj       ...      Xs] 

f    =    [/i     /j     ...     /m] 
0    =    [0     0     ...     0] 

We  call  numerical  methods  methods  because  in  general  they  are  not 
algorithms:  they  are  procedures  which  terminate  for  certain  inputs  with 
highly  specialized  mathematical  properties,  and  which,  when  they  do  ter- 
minate, terminate  in  time  only  loosely  boundable  by  the  structure  of  the 
input.  Termination  in  the  numerical  method  world  is  called  convergence. 

Conte  and  de  Boor,  in  their  chapter  on  systems  of  equations  and  un- 
constrained optimization  [7,  chapter  5],  describe  the  numerical  methods  we 
see  employed  in  the  systems  we  discuss  in  this  paper,  namely  the  methods 
of 

•  Linearizing  constraints,  in  which  we  approximate  a  non-linear  func- 
tion by  the  tangent  to  the  curve  of  the  function  at  a  particular  point, 
that  is,  by  the  line  with  the  slope  of  the  derivative  of  the  fimction 
evaluated  at  a  point  on  the  curve  and  passing  through  that  point. 


• 


Minimizing  least  squares,  in  which  the  norm  defined  by  the  square  of 
the  errors  of  successive  approximations  to  the  partial  derivatives  of  a 
function  is  minimized. 

Relaxation,  in  which  the  values  of  variables  are  incrementally  modified 
and  the  system  error  evaluated  to  determine  whether  an  improvement 
has  been  obtained. 

Newton- Raphson  iteration,  or  the  search  for  roots  by  refining  intervals 
indicated  by  the  sign-crossing  of  derivatives  of  the  polynomials  in  the 
set. 
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3.1      Linearizing  constraints 

3.1.1     How  to  linearize  constraints 

Given  a  system  of  equations  f  (x)  =  0,  it  is  sometimes  desirable  to  obtain 
a  linear  approximation  to  the  system  at  a  point  Xo- 

This  may  be  accomplished  through  a  trivial  application  of  Taylor's  The- 
orem from  calculus,  namely  that  for  a  continuous  function  f{x), 

/(x)«/{a)  +  /'(a)(x-a) 

Hence  a  linear  approximation  to  a  system  f  (x)  =^0  about  a  point  Xo  is 
simply 

f(x)«f(xo)  +  f'(xo)(x-xo) 

where  f  is  the  Jacobian  matrix  (system  derivative),  defined  by 


f'(x) 


3.1.2      Constraint  linearization  in  SKETCHPAD 

Every  type  of  constraint  relation  in  SKETCHPAD  is  represented  by  a  subrou- 
tine which  will  compute,  for  a  given  initial  set  of  values  for  the  constrained 
variables,  the  error  E  ,  that  is,  the  difference  between  the  current  value  of 
the  system  and  the  desired  value  of  the  system,  as  specified  by  the  con- 
straint. For  example,  a  line  is  represented  by  a  pair  of  points  (xi,yi)  and 
{x2,y2)-  The  error  for  the  constraint  that  the  line  be  horizontal  is  defined 
to  be  the  difference  of  the  y-coordinates  of  the  two  endpoints,  |yi  —  yjl- 

Because  the  constraints  are  user-definable  subroutines,  the  constraints 
do  not  have  to  be  linear.  "The  computation  of  the  error  may  be  non-linear 
or  time  dependent,  or  it  may  involve  parameters  not  a  part  of  the  drawing 
such  as  the  setting  of  toggle  switches,  etc."  [36,  pages  113-114]  Sutherland 
explains  the  art  of  defining  errors  on  constraints  ais  follows: 

In  order  to  make  the  constraints  work  well  together,  it  is  nec- 
essary that  they  be  balanced,  that  is  that  the  partial  derivative 
of  error  with  respect  to  displacement  be  nearly  equal  for  all 
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constraint  types.  I  have  arbitrarily  tried  to  make  the  error  sub- 
routines compute  an  error  about  proportional  to  the  distance  by 
which  a  variable  is  removed  from  its  proper  position.  In  other 
words,  many  of  the  existing  constraint  computation  subroutines 
maie  the  partial  derivative  about  unity. 

So  for  each  constraint  c,  in  the  system,  imagine  that  we  have  a  possibly 
non-linear  error  function  subroutine  /,(x)  which  gives  the  "displacement 
from  proper  position"  of  the  variables  related  by  c,.  Let  Xq  be  the  initial 
value  of  the  system  whose  current  state  is  x.  Let  E  be  the  linear  approxi- 
mation to  the  error  in  the  system  in  state  x,  so  that 

E    «    f(x) 
E,    «    /,(x) 

and 

Eio  «  /,(xo) 

is  the  error  in  the  initial  state  of  the  system.    Then  Sutherland  defines  a 
numerical  approximation  to  the  linear  approximation  of  /i(x)  as 

Ei  =  £.0  +  E  ^(=^.  -  ^O 


where  we  define 
and  note  that 

and 


AEi  =  Ei-  Ei, 
AjB,       dEi 


Axi        dxi 


Xo 


/i(x)  «  /i(Xo)  +  V/,|^  •  (x  -  Xo) 

The  general  idea  here  is  that  we  want  to  minimize  the  error  E,  that  is, 
we  are  trying  to  drive  the  expression 

Eo  +  VE-{x-  Xo) 

to  zero.  Given  that  we  now  have  a  linear  approximation  to  the  system,  the 
technique  of  minimizing  least  squares  which  we  discuss  next  will  help  us  do 
so. 
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3.2     Minimizing  least  squares 

The  linear  approximation  to  a  nonlinear  system  given  in  the  previous  sec- 
tion will  not  be  square  if  the  original  system  was  not  square.  In  this  case, 
there  may  be  no  exact  solution  as  the  resulting  linear  system  may  be  under- 
constrained  or  inconsistently  overconstrained.  To  get  around  this  problem, 
Sutherland  applies  the  least  mean  squares  approximation  method  to  in- 
erpolate  the  inconsistent  constraints  to  yield  a  compromise  square  linear 
system  which  may  then  be  solved  via  Gaussian  elimination.^ 

Suppose  we  have  approximated  a  nonsquare  nonlinear  system  f  (x)  =  0 
to  a  nonsquare  linear  system  fz,(x)  —  0  where 

h    =     [Ili     /lj     •••     Ilm] 

N 

Equivalently,  we  may  say  that 

h{X)  =  AX-C 

using  a  matrix  representation. 

Ideally,  /i,.(x)  =  0,  for  1  <  t  <  M .  The  error  in  /i,,,  jEx,.,  is  simply  the 
extent  to  which  /l,(x)  does  not  equal  zero: 

El,  =  /l.(x) 

Then  the  2-norm  of  Ex,  =  ^lIx),  El,,  is  a  measure  of  the  error  in  fL(x) 
in  which  we  square  the  error  of  each  component  of  ff,: 

When  we  minimize  the  2-norm  E^,  of  f^,  we  have  also  minimized  f^,. 
To  minimize  Ei,,  it  is  heuristically  sufficient  to  solve  for  the  case  in  which 


*The  following  exposition  is  adapted  from  Appendix  F  of  Sutherland's  thesis  |36],  with 
Bome  changes  in  notation. 
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the  gradient  of  El,  Js  0.  That  is,  we  set 


Since 


ox  ox 

(^  d         \ 


Now 


M 


1=1  " 


M 

E 

t=i 


^r 


J2  0.ijXj  —  Ci 


But 


5xt 


N 


So 


VE, 


l<Jk<A? 


Setting  VEl,  =  0,  we  want  to  find  x  such  that 


M 


N 


(N    M  M 

I31I(aita,y  -  Xj)  -  X)«.t<^i 
3=1 .=1  .=1 

In  matrix  notation,  this  may  be  expressed  as 

{A^A)X  =  A^C 


l<k<N 


l<k<N 
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Now,  A'^A  is  square  of  order  N,  hence  it  remains  to  solve  the  square 
system 

{A'^A)X-  A'^C  =  0 

for  X  via  Gaiissian  elimination. 

3.3      Relaxation 

3.3.1     Defining  relaxation 

Relaxation  is  a  numericzLl  method  for  approximating  the  zeroes  of  a  possibly 
nonlinear  square  system  of  equations  f  (x)  =  0.  Given  aji  initial  solution 
x^°\  Conte  and  de  Boor  [7]  describe  relaxation  as  follows. 

Suppose  that  the  /,  of  f  can  be  ordered  such  that  it  is  possible  to  find 
a  solution 

(7,(xi,.  ..,i,_i,x,+i,. . .  ,xs)  =  Xi 

of  /,(x)  =  0  for  1  <  t  <  TV. 

Then  given  an  absolute  error  bound  e  on  x^''^  —  x^'\  and  an  initial 
solution  x(°),  the  relaxation  procedure  is: 

function  RelaxationSolve  (g.  x(°^  e) 
begin  i  :=  0 
repeat 

i  :=  i  +  1 

for  i  :=  1  to  N 

ifi<7,(xW)-x;')|<|xW-x;-'^i 

then  xj'-^)  :=  x^p 

until  |xW  -x('-^)|  <  e 
return  x^'^ 
end 

In  other  words,  we  keep  replacing  components  of  x  by  solutions  of  f  (x) 
for  the  other  components  of  x  evaluated  at  the  current  value  of  x,  when 
this  does  not  diverge,  until  the  total  improvement  is  less  than  the  absolute 
error  boimd.  Each  pass  which  replaces  all  components  is  called  a  sweep. 
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3.3.2     Relaxation  in  SKETCHPAD  and  Thinglab 

As  an  example  of  the  combined  application  of  the  above  methods,  when  the 
network  technique  of  propagation  of  degrees  of  freedom,  described  above, 
does  not  work,  SKETCHPAD  performs  constraint  satisfaction  on  the  system 
f  (x)  =  0  in  a  given  initial  state  Xo  by 

1.  Linearizing  f  to  obtain  fi- 

2.  Making  fx,  square  by  the  least  mean  squares  method,  to  obtain  fx,^. 
.3.  Employing  releixation  on  fis  to  obtain  a  final  satisfying  assignment 

Xj7. 

Thinglab  also  calls  on  relaxation  to  solve  networks  of  relationships  in- 
volving simple  arithmetic  operators  and  equality  [3].  Borning  cites  Suther- 
land as  his  guide  in  the  implementation  of  the  relaxation  method.  Doming 
uses  the  technique  uses  the  technique  of  trimming  the  constraint  graph 
to  be  relaxed  (basically  the  network  technique  of  propagation  of  degrees 
of  freedom)  to  reduce  the  size  of  the  systems  presented  to  the  relajcation 
solver. 

3.4      Newton-Raphson  iteration 

3.4.1      Details  of  the  method 

To  understand  this  method,  let  us  first  consider  the  scalar  polynomial  case. 
Lipson  provides  a  nice  intuitive  explanation  of  Newton-Raphson  iteration 
in  this  caise  [22,  pages  313-314]. 

Consider  the  scalar  polynomial  equation  /(i)  =  0,  and  an  approxima- 
tion to  a  zero  of  /,  x*.  The  tangent  line  through  the  point  (it, /(it))  will 
intersect  the  i-zocis  at  a  point  (it+i,0)  which  is  closer  than  it  to  a  zero 
of  /.  The  slope  of  this  tangent  line  is  /'(lo)-  The  equation  of  the  tangent 
line  is  y  =  /'(it)i  -I-  b.  But  the  tangent  passes  through  (it, /(it)),  so 
b  =  f{xk)  -  f'{xk)xk-  Solving 

y  =  f'{x,)x  +  f{x,)-f'{x,)x. 
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for  the  i-axis  intersection  point  (xjk+i,0),  we  get 

0    =     f'{xt)xt+i  +  f[x,)-f'{xt)xt 


Xk+\      =      Xjk  - 


f'M 


The  last  equation  is  the  basic  Newton-Raphson  iteration  step.  Note 
that  an  important  alternative  expression  for  this  step  is 

m  +  I^Ax  =  0 

We  may  extend  Newton-Raphson  iteration  to  systems  of  equations.  The 
extended  algorithm  applies  to  smooth  vector  functions.  The  concept  of 
smoothness  is  famihar  to  differentiaJ  geometers.  We  repeat  the  definition 
here  for  the  benefit  of  non-specialists: 

A  function  /  is  eontinuously  diffcrenttable  on  an  interval  [a,  6]  when  its 
derivative  exists  on  (a,  6)  and  is  continuoxis.  A  continuously  differentiable 
function  is  called  C^.  Letting  /^"^  denote  the  n**'  derivative  of  /,  if  it  is 
the  caise  that  for  0  <  i  <  n,  /^'^  is  continuously  differentiable  and  f^"^  is 
continuous,  then  we  say  that  /  is  C".  If  /  is  C*  for  all  integers  k,  we  say 
/  is  C°°,  or  smooth  [10]. 

So  let  e  be  an  error  bound,  and  x^°^  a  close  enough  guess  to  the  zero  of  f, 
a  smooth  vector  function  such  that  f  is  continuous  and  f'{C)  is  an  invertible 
matrix  for  any  constant  vector  ^.  Then  Newton-Raphson  iteration  method, 
defined  a^  follows,  converges  on  a  zero  of  f: 

function  NewtonSolve  (f.  Xo.  e) 
begin  i  :=  0 
repeat 

x(.+i)_xW-r(x('))"'f(x(')) 

until    |x('-i)  -x^l  <  e 
e  :=  x(') 
return  ^ 
end 
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Note  that  an  alternative  expression  for  the  iteration  step,  corresponding 
to  the  alternative  expression  given  for  the  scalar  case,  is 

/(.,y,z,...)  +  f{A.+  |^Ay+|(A.  + 

Forming  this  expression  with  particular  values  of  x^'^  and  x^'"^^  where 
X  =  [ X  y  z  . . .],  and  then  applying  Gaussian  ehmination,  is  equivalent 
to  computing  the  expression  x^'^  —  f  (x^'^)     f(x(')). 

3.4.2     A  way  of  programming  Newton  iteration 

Wilkes  [39]  has  noted  that  Newton-Raphson  iteration  may  be  performed 
by  calling  at  compile-time 

•  An  underlying  algebraic  manipulation  system  for  performing  formal 
differentiation  (with  algorithms  such  as  those  explored  in  early  work 
by  McCarthy  [24]),  and  by  calling  at  rim-time 

•  A  linear  algebra  package,  for  numerically  solving  systems  of  linear 
equations. 

There  are  four  steps  to  the  process: 

•  Declare  names  of  constrained  variables  with  initial  approximations. 
These  are  called  "Newton  variables". 

•  By  formal  differentiation,  convert  statements  of  the  form  /(i,y)  =  0 
into  statements  of  the  form 

/(.,«) +|^A.+  |^A„  =  0 

•  For  a  block  with  declared  Newton  variables,  solve  the  linearized  sys- 
tem corresponding  to  the  variables,  \ising  increments  A3:,Aj/  and 
modify  i,  y  and  repeat  this  step  imtil  Ai,  Ay  are  less  than  some 
tolerance. 
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3.4.3     Constraint  systems  which  use  Newton-Raphson  iteration 

The  Juno  graphical  editor  uses  Newton-Raphson  iteration  to  solve  con- 
straints on  points  in  the  plane  arising  for  the  relationships  of  pairs  of  points 
being  parallel,  being  equal  in  length,  being  horizontal  and  being  vertical 
[27].  The  graphical  language  construct  for  defining  constrained  point  vari- 
ables requires  initial  solutions  for  their  location,  in  order  to  disambiguate 
solutions  in  nonlinear  cases.  Hence  this  language  is  similar  to  Wilke's  pro- 
posal, described  above.  Nelson,  like  Borning  in  the  case  of  relaxation,  uses 
the  technique  of  trimming  the  constraint  graph  to  be  relaxed  to  reduce  the 
size  of  the  systems  presented  to  the  Newton-Raphson  solver. 
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We  discuss  two  kinds  of  method  in  this  section: 

•  Relatively  ad  hoc  variations  on  Gaussian  elimination,  created  by  ex- 
perimental conaputer  scientists. 

•  Algebraic  techniques  for  the  order-preserving  approximate  solution 
of  systems  of  multivariate  polynomial  equations,  created^  by  mathe- 
maticians and  theoretical  computer  scientists. 

4.1      Variations  on  Gaussian  elimination 

Several  computer  scientists  have  embedded  implementations  variations  of 
Gaussian  in  their  constraint  satisfaction  systems.  Here  we  briefly  discuss 
three  similar  efforts. 

4.1.1      The  LELAND  Elimination  Algorithm 

The  general  idea  Van  Wyk  is  optimistic  about  the  prospects  for  a  simple 
elaboration  on  Gaussian  elimination.  In  his  thesis  [38,  page  35],  he  asks 

We  would  like  to  design  a  figure  l&nguage  with  a  few  powerful  primi- 
tives that  can  be  implemented  on  machines  of  moderate  size.  Do  we 
really  want  a  huge  symboUc  algebra  i>ackage  as  part  of  this  system? 
If  you  said  "no",  read  on  . . . 

He  chose  to  develop  a  compact  application-specific  solver.  This  solver, 
an  elimination  algorithm,  is  discussed  further  in  a  1984  paper  with  Emanuel 
Derman  at  AT&T  Bell  Labs  on  the  applications  of  this  solver  to  a  financial 
modelling  tool  [9].  A  rough  paraphrcise  of  the  algorithm  is: 

1.  Order  the  input  set  of  equations  ba^ed  on  their  arithmetic  complexity. 

2.  Perform  backward  substitutions  in  this  order  from  simplest  to  most 
complex. 


^Or  discovered,  depending  on  whether  one's  self-image  is  that  of  a  scientist  or  an  engineer. 
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3.  Simplify  the  resulting  set,  eliminating  tautological  equations. 

4.  Iterate  until  no  further  reduction  occurs. 

That  is,  do  roughly  what  you  would  do  with  paper  and  pencil,  taking 
some  care  to  make  sure  that  operation  on  a  polynomials  are  linear  in  the 
size  of  the  polynomial. 

They  find  that  their  solver  can  often  handle  simple  non-linear  equation 
systems  which  occur  in  practice,  for  exaimple  systems  which  contain  nonlin- 
ear equations  but  which  are  sufficiently  overconst rained  that  the  nonlinear 
equations  may  be  reduced  to  linear  ones  by  earlier  substitutions. 

Technical  details  Leler's  thesis  contains  a  nice  summary  of  the  details 
[21,  pages  33-37].  Derman  and  van  Wyk  present  the  precise  details  of  the 
algorithm  [9]. 

Imagine  an  input  set  of  multivariate  polynomials  equated  to  zero. 

Let  the  variables  be  totally  ordered,  for  example,  ordered  lexicograph- 
ically, into  a  sequence  Xi,  Xj, . . . ,  x^y^.  Define  am  ordered  linear  eombination 
to  be  a  normal  form  of  a  linear  equation,  in  which  we  cast  it  into  the  form 

N 

»=o 

Note  that  arithmetic  can  be  performed  on  ordered  linear  combinations 
as  follows: 

•  Multiplication  or  division  by  a  constant.  Simply  multiply  or  divide 
each  of  the  a,  by  the  constant. 

•  Addition  or  subtraction  of  two  ordered  linear  combinations.  Since 
the  terms  are  ordered,  simply  perform  an  "additive  merge"  of  the 
cofficients,  that  is,  just  add  the  coefficients  pairwise  to  obtain  the 
coefficients  of  the  result.  Similarly  for  subtraction. 

The  purpose  of  the  above  representation  is  to  be  able  to  do  solve  linear 
equations  in  linear  time  as  much  as  possible.* 


Mishra  notes  that  an  earlier  paper  due  to  Wu  may  also  have  ideas  on  this  topic. 
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Imagine  that  each  eqiiation  is  stored  in  a  parse  tree,  in  which  the  interior 
nodes  are  arithmetic  operators  such  as  +,  —  and  x,  and  the  leaves  are 
variables  and  integer  constants.  Linear  subexpressions  of  this  tree  can 
be  transformed  into  ordered  linear  combination  form  in  a  fairly  obvious 
fashion,  in  linear  time.  Call  this  the  normal  form  of  the  equation. 

The  main  algorithm  processes  the  set  of  the  normal  forms  of  the  equa- 
tions by  the  following  rules: 

•  If  the  equation  evaJuates  to  zero  directly  as  a  result  of  coefficient 
simplification,  discard  it. 

-  •  If  the  equation  evaluates  to  a  non-zero  constant,  then  fail:  the  system 
is  inconsistent. 

•  If  the  equation  is  imivariate,  say  6  +  ai  =  0,  then  substitute  —b/a  for 
X  in  all  the  remain  equations,  and  delete  this  equation  from  the  set, 
storing  the  solution  for  z  in  a  table. 

•  If  the  equation  is  multivariate  but  contains  a  term  which  is  univariate, 
solve  for  that  term,  and  substitute  for  that  term  in  all  other  equations. 
Store  the  expression  e  for  the  substituted  variable  z  in  a  table;  if  the 
remaining  equations  are  then  solved,  substitute  the  solutions  back 
into  e  to  solve  for  z. 

4.1.2     Algebraic  transformationB  in  Magritte 

Gosling's  Magritte  editor  employs  a  simplifier  of  the  term  rewriting  sys- 
tem variety  to  simplify  algebraic  expressions  arising  from  constraint  graphs 
with  arithmetic  relations  and  numerical  variables.  Like  van  Wyk  (see  Sec- 
tion 4.1.1),  he  maintains  equations  in  a  cjinonical  form;  in  Gosling's  case, 
the  sum  of  products  of  factors. 

This  simplifier  operates  with  17  rewrite  rules  on  expressions  involving 
-H,  — ,  X,  -r,  exponentiaion  and  equality,  which  he  claims  are  sufficient  for 
almost  all  of  the  practical  constraint-resolution  problems  he  encounters. 
For  the  record,  they  are 

-c    =>     -U 
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«i  -  ej    =>    «!  +  (-lej) 

Ci  -r  Cj      =>      ^i^j 
61(62  +  63)      =>      6162  +  CiCs 

el  =>  e 

e+0  =>  6 

kie  +  k2e  =>  (^1  +  ^2)6 

e'  =>  6 

c°  =>  1 

(6162)'^      =>      6^65' 

e\    =>    6161 
el    =>    616161 

61  =  62      =>      61  —  Cj  =  0 
61  A  Ci     =>     61 
61  V  61     =>     61 

In  an  appendix  on  his  simplifier,  he  exclaims  ° 

It  IB  very  simple,  certainly  much  simpler  than  MaCSYMA.  But  that 
is  precisely  the  point:  powerful  results  can  be  had  by  applying  even  a 
very  simple  amoimt  of  algebraic  knowledge. 

4.1.3     Leler's  augmented  term  rewriting  system 

William  Leler  defines  an  augmented  term  rewriting  system  language,  BeRT- 
RAND,  that  is  a  production-system  like  language  plus  single  assignment  to 
variables  and  types  with  single  supertypes  [21] 

This  language  is  Turing  equivalent.  Under  certain  restrictions  on  the 
niles  given  in  the  thesis  [pages  43-44],  fast  pattern  matching  algorithms 
may  be  used  for  rule  rewriting. 


"This  remark  b  criticiied  by  Yap,  who  informed  the  author  in  conversation  in  October 
that,  while  the  remark  hae  its  place,  Gosling's  thesis  was  too  impressionistic  from  a 
technical  point  of  view,  to  give  Gosling  a  basis  to  make  the  comment. 
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It  provides  a  convenient  language  for  implementing  algebraic  algorithms, 
with  the  caveat  that  the  user  has  much  less  command  over  the  flow  of  con- 
trol than  with  a  conventional  language.  This,  however,  is  claimed  as  a 
benefit,  as  the  term  rewriting  system  mechanism  control  structure  often 
results  in  the  ability  to  automatically  invert  functions  without  reprogram- 
ming. 

The  paradigm  of  augmented  term  rewriting  systems  is  strikingly  similar 
to  context  senstive  grammars  and  their  derivations,  and  also  quite  similar 
to  the  syntactic  conversion  rules  of  Church's  A-calculus. 

For  further  discussion  of  the  hnguistic  sispects  of  BERTRAND,  see  Sec- 
tion 6.8. 

4.2      Advanced  algebraic  methods 

There  are  several  techniques  from  the  theory  of  algorithmic  algebraic  ge- 
ometry which  may  be  of  use  to  us. 

Given  that  our  general  problem  is  the  solution  of  systems  of  polynomial 
equations,  we  start  with  a  few  results  on  the  complexity  of  subcases  of  this 
problem.  This  is  not  a  comprehensive  survey,  but  rather  we  plaice  these 
observations  here  as  a  reminder  of  the  potential  difficulties  and  of  the  need 
to  restrict  the  langiiage  of  Section  1.2  to  the  least  powerful  subset  consistent 
with  interesting  research. 

We  discuss  a  particular  combined  application  of  two  methods,  Grobner 
bzises  and  algebradc  number  arithmetic,  which  will  allow  us,  in  theory,  to 
efficiently  compute  zeroes  of  systems  of  polynomial  equations  of  interest 
to  us,  to  an  order-preserving  degree  of  accuracy.  Note  that  Grobner  bas- 
es have  been  more  intensively  investigated  experimentally;  we  even  give  a 
short  example  of  a  running  program! 

4.2.1     van  Wyk's  observations  on  algebraic  complexity 

The  following  observations  are  a  summary  of  a  summary  from  van  Wyk's 
thesis  [38,  Appendix  l]. 

Let  F  be  a  set  of  variables,  £^  be  a  set  of  simultaneous  equations,  and  a, 
6,  and  c  be  variables  in  V  or  rational  numbers.  For  the  first  three  results, 
the  problem  is  to  find  a  satisfying  assignment  to  E: 
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•  SimuHantous  maximization  (NP  Complete).  E  consists  of  equations 
of  the  form  max(a,  6)  =  c. 

•  Simtiltaneous  absolute  values  (NP  Complete).  E  consists  of  linear 
combinations  and  absolute  values  of  lineair  combinations. 

•  Simultaneous  quadratics  (NP  Complete).  E  consists  of  equations  of 
the  form  {a  +  b){c  +  d)  =  0. 

•  Sin-abs-polynomiaJ  zero  equivalence  (Recursively  undecidable).  Given 
an  expression  e  involving  r,  +,  — ,  sin,  abs  and  composition,  does 
e  =  0? 

4.2.2      Grobner  bases  and  ideal-theoretic  methods 

Grobner  bases  are  a  bases  for  poljmomial  ideals,  with  computational  prop- 
erties that  suit  them  to  a  host  of  decision  problems,  and  an  elegant  un- 
derlying theory  [26].  In  particular,  there  is  a  technique  involving  Grobner 
bases  which  leads  us  towards  the  solution  of  systems  of  polynomial  equa- 
tions [5].  Let  us  explore  the  Grobner  basis  technique  only  as  it  applies  to 
this  problem. 

Definition  of  Grobner  bases.  Let  F  be  a  system  of  algebraic  equations, 
that  is,  a  finite  subset  of  R  =  /r[ii,i2, . . . ,  i„]  where  if  is  a  field.  We  may 
view  any  system  of  equations  as  the  basis  of  an  ideal.  For  example,  R 
could  be  Q[x,y,z]  in  the  case  of  systems  of  equations  in  three  variables 
with  rational  cofficients,  or  Q{a,b)\x,y],  in  the  case  of  linear  equations 
with  coefficient  expressions  in  the  "symbolic  constants"  a  and  6,  and  let  g, 
hi  and  /ij  denote  polynomitils  in  F. 

The  set  of  equations  generated  by  the  ideal  will  have  the  same  set  of  ze- 
roes as  the  original  system.  Several  questions  may  be  easily  answered,  such 
as  the  existence  of  a  finite  or  infinite  number  of  solutions,  if  we  can  trans- 
form our  ideal  into  a  Grobner  basis.  Further  speedup  and  space  economy 
may  result  if  we  compute  a  reduced  Grobner  bzLsis  [26]. 

A  monomial  is  a  term  of  the  form  ap  where  a  is  a  constant  and  p 
is  a  product  of  powers  of  Xx,i2, . .  .,x„.    Let  Hmono(/)  denote  the  head 
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monomial  of  a  polynomial  of  /  under  some  fixed  total  ordering  (such  as 
lexicographic)  of  products  of  powers  of  variables. 

A  reduction  of  g  with  respect  to  F  is  a  reduction  of  g  with  f  ^  F  defined 
to  be  the  result  h  given  by 

h  =  g-cf 

where 

m 
c  = 


Hmono(/) 

and  m  is  some  monomial  in  g. 

A  normal  form  of  g  is  defined  to  be  the  transitive  closure  of  the  reduction 
of  g  with  respect  to  polynomials  in  F. 

(The  above  details  follow  Mishra  and  Yap's  notation  [26]). 

Then  Buchberger  defines  a  F  to  be  a  Grobner  basis  when  for  all  g,  hi,  hi  G 
F,  if  hi  and  h^  are  normal  forms  of  g  mod  F  then  hi  —  h^. 

A  sketch  and  example  of  the  intended  algorithm.  A  Grobner  basis 
G  derived  from  an  initial  basis  (F)  has  the  property  that,  if  1  is  in  the 
basis,  then  the  system  is  inconsistent  [5,  Method  6.8]. 

There  is  a  fast  method  for  determining  whether  G  has  infinitely  many 
solutions  [5,  Method  6.9]. 

G  has  the  property  that,  if  ij,  Xj, . . .  i„  are  the  variables  in  F,  then  one 
polynomial  in  G  will  be  univariate.  Without  loss  of  generality,  asstmie  that 
polynomial  is  univerate  in  ii.  There  will  be  n  —  1  remaining  polynomials, 
each  with  one  term  x,,  and  the  remaining  terms  in  Xi  [5,  Methods  6.10, 
6.12]. 

We  may  now  observe  that  the  problem  reduces  to 

1.  The  problem  of  finding  zeroes  for  single  univariate  polynomial,  fol- 
lowed by 

2.  The  performance  of  a  series  of  substitutions  into  the  remaining  poly- 
nomials, which  are  expressed  as  the  diff^erence  of  one  of  the  oaultiple 
variables  and  a  polynomial  in  the  first  variable. 

For  this  last  step,  we  would  like  to  be  able  to  compute  an  approximate 
representation  of  the  zeroes  of  the  system,  which  is  simuitaneously  order 
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preserving  in  the  sense  that,  if  there  are  a  sum  of  m  zeroes  for  all  the  equa- 
tions in  G,  then  for  no  two  zeroes  will  it  be  the  case  that  our  approximate 
representation  allows  us  to  falsely  conclude  their  numerical  order.  If  we  are 
using  algebraic  numbers  (see  Section  4.2.3),  this  means  that  we  expect  the 
defining  intervals  to  be  pairwise  disjoint. 

How  to  compute  a  Grobner  basis.  Mishra  and  Yap  [26]  discuss  the 
complexity  of  the  algorithm  for  computing  the  Grobner  basis  of  a  system 
F.  Of  course,  Buchberger  also  discusses  the  matter,  and  gives  several  im- 
provements over  the  basic  method  [5,  Methods  6.10,  6.12]. 

Let 

m  m 

^^^'^'  "  Hmono(/)-^  "  Hmono(ff)^ 
Then  Mishra  and  Yap's  presentation  of  the  basic  method  is  as  follows: 

Algorithm  GrobnerBasis  {F) 

Input  A  finite  set  of  polynomials  F. 

Output  A  Grobner  basis  G  for  F. 

Method 
begin 

G.=  F 

B~{{f,9}:f,9eGJj^g} 
whWeB  ^9  6o 

Choose  {f,g}  from  B;  B  :=  B  -  {f,g} 
h:=S{f,g) 

h'  :=  the  normzd  form  of  h  with  respect  to  G 
if/i'  /  0  then 
B:=Bu{{f',h'}:f'eG} 
G:=Gu{/i'} 
end  if 
end  while 
return  G 
end 


4.2     Advanced  algebraic  methods  45 

An  application  from  SCRATCHPAD- II.     As  a  demonstration  of  both 

•  The  power  of  an  implementation  of  the  Grobner  basis  method  to 
simplify  the  programming  involved  in  constraint  satisfaction,  and 

•  The  program  complexity  of  a  simple  constraint-satisfaction  problem 
(taking  a  dumb  approach  and  a  fancy  one), 

Consider  the  problem  of  constraining  the  location  of  a  point  to  lie  at  the 
intersection  of  two  lines.  Six  variables  are  involved  here:  the  slope  and  y- 
intercept  of  the  two  lines,  and  the  x  and  y  coordinates  of  the  point.  Given 
any  four  of  the  quantities,  one  can  form  a  system  of  two  linear  equations, 
the  solution  of  which  yields  the  other  two.  One  way  to  solve  this  problem 
in,  say,  the  language  PASCAL,  without  doing  algebraic  expression  manip>- 
ulation,  is  to  enumerate  the  f^j  possible  combinations  of  known  variables, 
and  then  for  each  of  these  fifteen  cases  write  a  piece  of  code  which  solves 
the  appropriate  linear  system.  So  a  procedure  in  PASCAL  for  a  "line-line 
intersection  constraint  node"  would  require  perhaps  100  Unes  of  code.  The 
following  code  of  Rudiger  Gebauer  [12]  shows  how  to  solve  this  problem 
with  his  Grobner  basis  package  for  the  SCRATCHPAD  II  computer  algebra 
system[l4].  First  we  defijie  a  prototype  linear  system: 

1.  pi  :  DMP[i,  y]  RF  RN  :=  Cii  +  Cjy  +  Cj 

2.  pj  :  DMP[i,  y]  RF  RN  :=  6ii  +  ftjy  +  63 

These  are  two  assignment  statements.  The  expressions  involving  the  DMP, 
RF  and  RN  are  types,  indicating  that  variables  pi  and  pj  are  distributed 
multivariate  polynomials  in  x  and  y  over  rational  fields  of  rational  numbers. 

3.  Pn  ••=  [Pl,P2] 

Here  we  are  just  forming  a  list  of  n  polynomials,  for  n  =  2.  Now  we  invoke 
the  Grobner  basis  package. 

4.  groebner(p„) 
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This  causes  the  system  to  print  the  following  expression  pair: 
t        — ojfts  +  0362  0163  —  as6i , 

F  +  — I r~'  y  +  ~~k rl 

These  are  zeroes  for  i  and  y,  that  is,  solutions. 

We  have  asked  for  the  result  in  terms  of  x  and  y.  To  get  the  result  for, 
say,  as  and  62,  presumably  we  just  list  03,62  in  place  of  x,y  in  the  type 
expressions  in  lines  1  and  2,  and  the  package  does  the  rest.  Bjised  on  the 
example  as  it  stands,  the  SCRATCHPAD  II  procedure  is  five  lines  of  code: 
the  above,  plus  a  procedure  header  which  takes  as  its  argument  the  pair  of 
variables  to  be  substituted  into  the  type  expressions  in  lines  1  and  2. 

4.2.3     Algebraic  numbers 

A  complex  number  ^  is  called  an  algebraic  number  if  it  is  a  zero  of  a  poly- 
nomial equation  /(i)  =  0  where  /  €  Q[x]  [29].  (We  are  only  interested  in 
real  algebraic  numbers.) 

The  utility  of  algebraic  numbers  is  that  we  can  perform  order-preserving 
arithmetic  on  approximations  to  them.  In  particular,  this  allows  us  to  do 
arithmetic  with  the  solutions  of  polynomials.  (For  an  application  relevant 
to  the  topic  of  this  paper,  see  Section  4.2.2.)  By  order  preserxnng  we  mean 
that  if  two  algebraic  numbers  ^1  and  ii  are  such  that  ^1  <  ^2>  then  it 
is  possible  to  refine  the  isolating  intervals  for  the  corresponding  interval 
representation  (see  below)  of  these  algebraic  nimabers,  say  ^j  and  ^j,  until 
it  is  the  case  that  ^j  <  ^j. 

The  arithmetic  algorithms  we  desire  on  approximate  algebraic  numbers 
are  just  comparison,  addition,  subtraction,  multiplication  and  division.  The 
algorithms  depend  on  certain  intermediate  results,  namely 


• 


A  method  of  isolating  all  roots  of  a  polynomial  baised  on  the  size  of 
the  coefficients. 

A  method  of  isolating  aJl  individual  roots  given  the  interval  of  all 
roots  of  a  polynomial.  We  discuss  the  method  of  Sturm  sequences 
[42],  which  has  an  0(n  log' n)  complexity,  where  n  =  deg[p).  We 
might  exclude  the  method  of  Collins  and  Loos  [41],  as  it  appears  to 
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have  a  much  higher  complexity,  namely  0{n^°  +  n7  log'  (f),  where  d  is 
the  length  in  bits  of  the  representation  of  p.  However,  good  average- 
case  performance  has  been  observed  for  the  latter  method  [43]. 

Loos  discusses  approximate  algebraic  number  arithmetic  in  the  context 
of  computing  in  algebraic  extension  fields  [23].  Schwartz  and  Sharir  also 
discuss  the  matter  in  an  appendix  to  a  paper  on  topological  computations 
in  robotics  [30,  Appendix  B]. 

Approximate  algebraic  numbers     Define  the  interval  representation  of 
algebraic  numbers  to  consist  of 

1.  A  defining  polynomial  of  the  form 

i=0 

where  the  a^  are  integers. 

2.  An  isolating  interval  consisting  of  pair  of  binary  rational  numbers 


W'2<^) 


where  a,  6,  c  and  d  are  integers. 

The  isolating  interval  is  capable  of  isolating  a  real  number  to  any  desired 
degree  of  precision,  due  to  the  fact  that  the  rational  numbers  are  dense,  that 
is,  that  between  any  two  reals  there  is  a  rational.  The  particular  rational 
number  representation  we  have  chosen  is  advantageous  in  the  sense  that 
halving  is  easy  and  results  in  a  rational  in  the  same  form. 

Essential  operations.     There  are  two  essential  operations  for  computing 
with  algebraic  numbers: 

•  Computing  the  bounding  interval  of  all  roots  of  a  polynomial. 

•  Isolating  the  individual  roots  within  the  bounding  interval. 
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Arithmetic  on  algebraic  numbers.     We  need  to  be  able  to  perform 
comparison  (a  direct  consequence  of  the  essential  operations  discussed  above), 
as  well  as  addition,  subtraction,  multiplication  and  division.    One  should 
consult  Loos'  paper  for  the  algorithmic  details. 
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5      Finite-precision  arithmetic  and  graphics 

For  maximum  efficiency  one  must  use  the  finite-precision  arithmetic  oper- 
ations available  on  a  typical  computer.  The  task  is  then  to  keep  a  record 
of  round-oflf  errors,  and  never  come  to  a  conclusion  which  presiunes  more 
information  than  is  actually  available.  This  approach  is  acceptable  in  prob- 
lems involving  polygonal  regions,  such  as  VLSI  layout,  but  does  not  seem 
acceptable  in  geometrical  computations  involving  points  with  no  area  and 
lines  with  no  width. 

Segal  and  Sequin  introduce  fornas  of  objects  with  minimum  feature  size 
and  face  thickness  of  three-dimensional  objects,  and  show  how  to  compute 
in  finite  precision  "with  topological  immunity"  on  objects  in  this  form  [31]. 
The  minimum  feature  size  depends  on  the  objects  dimensions  and  location 
in  space,  and  the  face  thickness  depends  on  "how  well  a  face's  vertices 
conform  to  its  computed  plane". 

Milenkovic  defines  two  techniques  [25]: 

•  Data  normalization,  which  transforms  a  geometric  structure  into  a 
structure  for  which  finite  precision  calculations  will  yield  only  topo- 
logically  correct  answers. 

•  The  hidden  variable  method  constructs  geometric  structures  which 
belong  to  an  infinite  precision  domain,  but  which  are  represented 
finitely. 

We  have  not  had  the  opportunity  to  review  a  recent  paper  by  Yao  and 
Greene  which  is  also  relevant  [40]. 
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6      Constraint  languages 

In  this  section  we  review  work  on  'programming  systems"  which  employ 
the  constraint  model  directly.  By  programming  systems  we  mean  either  pro- 
gramming languages  and  their  runtime  environments,  or  graphical  editors 
with  a  significant  embedded  linguistic  component. 

6.1      The  Sketchpad  graphical  editor 

Ivan  Sutherland  defined  an  interactive  graphical  editor  which  allowed  the 
definition  of  hierarchical,  structured  objects  with  simple  arithmetic  con- 
straints on  the  parts  of  the  objects.  His  SKETCHPAD  program  on  the  MIT 
Lincoln  Laboratories  TX-2  computer  in  1963  constituted  a  prototype  per- 
sonal computing  tool  which  was  about  15  years  ahead  of  its  time.  His 
use  of  constraints  was  novel,  and  this  use  prefigured  the  incorporation  of 
constraints  in  progranmaing  languages. 

Sketchpad's  graphical  programming  language  is  not  easy  to  describe 
syntactically,  a£  it  is  encoded  in  series  of  special  pushbuttons  and  light 
pen  motions  used  to  denote  particular  commands.  However,  it  is  easy  to 
analyze  the  overall  structure  of  the  language  for  which  the  SKETCHPAD 
graphical  editor  is  the  interpreter. 

Objects.     The  objects  manipulated  by  SKETCHPAD  include: 

•  Line  segments.  Denoted  by  two  endpoints  (u,v). 

•  Circular  arcs.    Denoted  by  a  center  point  u,  a  radius  r,  and  an  arc 
angle  9. 

•  Attachment  points. 

•  Digits. 

•  Text. 

•  Scalars. 

•  Dummy  variables. 
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Commazids.     Rotation,  positioning,  scaling,  recursive  merging. 

Constraints.     There  are  seventeen  geometric  constraints  types  in  SKETCH- 
PAD [36,  Appendix  A].  Omitting  some  obscure  ones,  these  are: 

•  Three  things  collinear 

•  Distance  from  first  to  second  equals  distance  from  first  to  third. 

•  Thing  is  erect  or  on  its  side. 

•  First  thing  directly  above,  below  or  beside  second  thing. 

•  One  thing  parallel  to  line  between  two  other  things. 

•  Distance  from  first  thing  to  second  is  1/3,  1/2,  1,  2,  3  times  distance 
from  third  to  fourth. 

•  First  thing  is  1/3,  1/2,  1,  2,  3  times  size  of  second  thing. 

•  Value  of  displayed  scalar  equals  distance  between  two  things  in  inches. 

•  Value  of  displayed  scalar  equals  size  of  thing  in  inches. 

•  First  thing  is  at  midpoint  of  other  two. 

•  Thing  is  1/32,  1/16,  1/8,  1/4,  1/2  or  1  inch  in  total  size. 

•  A  line  from  first  thing  to  second  thing  would  be  parallel  or  perpen- 
diciilar  to  a  line  from  third  thing  to  fourth  thing.  (The  lines  need  not 
exist.) 

Types.  Subpictures,  generic  structure,  definition  copying,  definitions 
(e.g.  horizontal  line,  dimension  line  with  arrowhead  and  number),  poly- 
morphic functions. 
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6.2      A  proposal  of  Wilkes 

In  1964  Maurice  V.  Wilkes  suggested  adding  equations  as  statements  to 
ALGOL-like  languages.  The  compiler  would  then  generate  code  for  a  New- 
ton iteration  solution  for  the  xmknown  variables.  He  proposed  to  add  two 
new  notations  to  ALGOL-Uke  languages: 

•  A  notation  for  declaring  variable  names  and  initial  approximations  to 
variable  values,  used  as  starting  points  for  Newton  iteration: 

c  ~  5 
X b/a 

•  Equations,  such  as 

25  X  cos(A)  +  45  X  cos(B)  =  60 
AxXt2+BxX+C=0 

Wilke's  suggestion  seems  original  and  qiiite  imaginative  for  the  time  (he 
was  programming  in  AUTOCODE  for  the  Edsac-2).  But  Prof.  James  Dem- 
mel  of  Courant  Institute  said  in  a  conversation  in  September  that  the  idea 
of  enhancing  an  ordinary  language  compiler  to  "automatically  program" 
loosely-specified  nimaerical  analysis  computations  has  been  developed  since 
then  under  such  names  as  "scientific  workbench" ,  for  example  at  the  Math- 
ematics Research  Center  at  the  University  of  Wisconsin.  The  approach  has 
fallen  into  disfavor,  because  anybody  who  really  wants  to  do  a  numerical 
analysis  computation  is  likely  to  want  to  have  exact  control  over  the  meth- 
ods chosen,  in  order  to  ensure  the  desired  convergence  level  for  a  particular 
application.  We  don't  know  enough  yet  to  design  a  language  with  enough 
application  information  in  it  to  tell  an  "expert  compiler"  what  it  needs  to 
know  to  select  among  the  many  numerical  methods.  It  is  possible  that,  as 
a  specialist  in  numerical  analysis,  he  is  speaking  as  one  who  wants  control 
over  the  finest  details.  For  the  applications  of  many  naive  users,  perhaps 
the  automatic  programming  approach  would  be  preferable.  For  example, 
this  paper  is  written  using  the  lATgpC  system,  which  makes  many  stylis- 
tic decisions  for  the  user  [18],  but  more  particular  user  may  prefer  naked 
1^[16]. 
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6.3      Constraint  languages  and  applications  at  MIT 

In  the  mid-to-late  1970's,  Gerald  Sussman,  Richard  M.  Stallman  and  Guy 
Louis  Steele,  Jr.  at  the  MIT  Artificial  Intelligence  Laboratory  began  a  line 
of  research  in  the  application  of  constraint  satisfaction  as  a  language  and 
computational  model  to  solve  problems  in  the  construction  of  computer- 
aided  design  systems  for  electrical  circuitry,  and  for  the  modelling  of  simple 
physical  systenis  such  as  bridges  imder  load  (Sutherland's  example)  and 
temperature  converters.  This  work  grew  out  of  similar  work  iising  a  logic 
programming  paradigm.  Stallman  and  Sussman  wrote  a  program  to  an- 
alyze circuits  in  1975  [35].  Steele  and  Sussman  wrote  a  joint  paper  on  a 
progranmiing  language  based  on  the  formation  and  solution  of  constraint 
networks  [32],  which  culminated  two  years  later  in  Steele's  thesis  on  this 
topic  [33]. 

Steele's  thesis  investigated  the  problem  of  retraction  in  a  constraint 
graph  which  is  being  edited  and  iteratively  resolved,  and  the  problem  of 
explanation,  or  tracing  of  the  computation  process  by  which  a  satisfying 
assignment  is  found  for  a  constraint  graph.  To  aid  the  human  editing 
the  constraint  graph,  he  introduced  a  lattice  of  assertable  values,  which 
under  various  circumstances  took  precedence  over  each  other  depending  on 
what  was  known  at  the  time.  The  basic  constituents  of  the  lattice,  from 
weakest  to  strongest  precedence,  are  defaults,  assumptions,  constants  and 
solutions.  The  latter  type  of  values  are  those  arrived  at  by  inspection  of 
constraint  relations,  constants,  assumptions  and  defaults.  The  history  of 
decisions  in  the  satisfying  jLSsignment  computation  algorithm  is  kept  for 
each  computed  value.  This  history  is  used  to  "undo"  or  incrementaJly 
recompute  the  satisfying  eissignment  when  assumptions  and  defaults  are 
retracted. 

Some  students  of  MIT  Professor  Richard  Zippel,  Farbus  and  Williams, 
have  embedded  constraint-satisfaction  techniques  in  a  VLSI  design  system 
[44].  Zippel's  Schema  design  system  has  embedded  in  it  a  Grobner  bas- 
is package  for  analyzing  the  properties  of  VLSI  circuits.  This  is  a  large, 
practical  VLSI  design  system  which  Zippel  claims  is  a  rival  to  the  MAGIC 
routing  editor  produced  at  U.C.  Berkeley.  We  may  taie  this  to  mean  that  it 
has  good  graphics,  solves  many  routine  placement  and  design  rule  problems, 
and  runs  fast. 
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6.4      The  THINGLAB  graphical  simulation  language 

Also  in  the  late  1970s,  Alan  Homing,  working  at  Stanford  University  and 
Xerox  Palo  Alto  Research  Center,  extended  the  SMALLTALK  language  to 
create  the  THINGLAB  constraint-oriented  simulation  program  [3].  He  did 
this  work  as  a  member  of  the  Xerox  PARC  Learning  Research  Group.  His 
aim  was  to  construct  a  language  and  programniing  environment  with  which 
it  would  be  particularly  easy  for  a  naive  user  (such  as  a  child)  to  statically 
describe  a  situation  and  then  explore  that  situation  through  simulation. 

Borning's  system  allowed  the  definition  of  new  types  of  object  and 
constraints  on  their  components.  He  developed  a  method  of  compiling 
constraint-satisfaction  code  which  did  local  propagation  of  values  and  com- 
putation of  unknown  values  based  on  simple  inversion  of  expressions  for 
known  values,  and  which  as  a  last  resort  invoked  Newton  iteration.  The 
details  of  his  constraint  satisfaction  methods  are  quite  similar  to  those  of 
Steele  (his  thesis  predates  Steele's  by  one  year). 

As  programming  language  design  Borning's  work  seems  superior  to 
Steele's,  building  as  it  did  on  the  rich  environment  of  Smalltalk,  with 
its  strong  graphical  orientation  and  attractive  clziss-oriented  type  system. 
Steele's  work  was  implemented  in  Lisp  Machine  Lisp.  The  Lisp  Machine  at 
that  point  did  not  have  an  attractive  system  for  representing  \iser-defined 
types. 

Borning's  system  is  embedded  in  the  linguistically  distinctive  SMALL- 
TALK 76  language,  which  allows  the  user  to  place  types  of  objects,  together 
with  procedures  for  operating  on  those  types,  in  a  hierarchy  of  subtypes 
(caJled  classes),  which  is  such  that  if  a  procedure  is  not  defined  for  a  sub- 
type but  is  defined  for  a  supertypxe,  then  the  supertype's  definition  may  be 
invoked.  A  procedure  definition  for  a  subtype  may  invoke  an  identically- 
named  procedure  in  the  supertype.  Procedures  in  SMALLTALK  are  called 
methods.  In  addition,  Borning's  version  of  SMALLTALK  was  extended  to 
allow  a  type  to  have  several  supertypes.  This  language  model,  generally 
referred  to  as  object-oriented  programming  [34],  is  not  of  major  concern  to 
us  here,  because  descriptive  paradigms  will  not  affect  the  complexity  of 
constraint-satisfaction  algorithms. 

The  following  example  shows  the  definition  of  a  line  with  a  point  con- 
strained to  lie  at  its  middle.   Note  that  the  constraint  is  expressed  as  an 
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equation  together  with  solutions  for  each  of  the  variables  in  the  equa- 
tion in  terms  of  the  other  variables  in  the  equation.  Also,  assume  that 
a  point  in  THINGLAB  is  represented  as  a  pair  of  integers,  and  a  line  as 
a  pair  of  points,  and  that  vector  arithmetic  on  points  is  implicit  (that  is, 
(0,6)  +  {c,d)  =  {a  +  c,b-\-d)): 

Class  MidPointLine 
Superclasses 

GeometricObject 
Part  Descriptions 
line:  a  Line 
midpoint:  a  Point 
Constraints 

(linepointl  +  line  point2)  /  2  =  midpoint 

midpoint  •<—  (line  pointl  +  line  point2)  /  2 

line  pointl  <—  line  point2  +  ((midpoint-line  point2)*2) 

line  point2  *—  line  pointl  +  ((midpoint-line  poinl)*2) 

Another  example  is  the  definition  of  a  triangle  and  an  Isosceles  trian- 
gle. This  example  introduces  the  concept  of  a  merge,  which  is  first  seen  in 
Sketchpad,  namely,  a  special  kind  of  constraint  which  asserts  the  equal- 
ity of  two  parts  (or  equality  of  variables,  or  equality  of  some  subset  of 
components  of  structured  variables): 

Class  Triangle 

Superclasses 

GeometricObject 
Part  Descriptions 

sidel:  a  Line 

side2:  a  Line 

side3:  a  Line 
Merges 

sidel  pointl  =  side3  pointl 

sidel  point2  =  side2  pointl 

side2  point2  =  side3  point2 

Class  Isosceles  Triangle 
Superclasses 
Triangle 
Constraints 
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sidel  length  =  side2  length 

6.5  The  MAGRITTE  graphical  editor 

Further  refining  the  line  of  applying  constraint  satisfaction  to  the  problem  of 
coniputing  of  satisfying  aissignments  to  constraint  graphs  with  simple  two- 
dimensional  representations,  James  Gosling  at  Carnegie-Mellon  University 
[13]  developed  MAGRITTE,  a  SKETCHPAD-like  editor  for  objects  in  the 
plane.  See  Section  4.1.2  for  a  description  of  the  algebraic  manipulation 
engine  within  MAGRITTE. 

6.6  The  Ideal  typesetting  language 

Christopher  John  Van  Wyk,  in  a  1980  Stanford  thesis  proposed  the  novel 
equational  language  IDEAL  for  typesetting  graphics.  IDEAL  provides  a 
block-structured  method  of  declaring  systems  of  equations  constraining 
variables.  He  introduced  a  somewhat  novel  block-structured  constraint 
language  which  is  reminiscent  of  Wilke's  proposal  (see  Section  6.2),  and 
which  is  also  a  precursor  of  the  JUNO  language  (see  the  following  section). 
The  language  consists  of  procedures  in  which  the  parameters  are  declared 
variables  for  which  the  user  may  or  may  not  supply  expressions.  The  proce- 
dures are  called  boxes.  A  box  call  consists  of  put  followed  by  the  box  name 
followed  by  defining  equations  for  the  parameters  the  user  wishes  to  specify 
in  brakes.  A  definition  for  a  box,  in  this  case,  for  a  circle  with  three  refer- 
ence points  for  attaching  lines  at  equidistant  points  on  its  circumference,  is 
as  follows. 

box  treenode  { 

var  rad.  cen.  height, 

hook.  Ichild.  rchild; 
rad  =  height  /  2; 
hook  =  cen  +  (O.rad): 
Ichild  =  cis{3*pi/4)[cen.hook]: 
rchild  =  cis(-3*pi/4)[cen.hook]: 
put  cirde{  radius  =  rad;  center  =  cen;  };  } 

The  idea  is  that  if  the  user  specifies  a  defining  expression  for  hook,  then 
this  is  equated  to  cen  -f-  (O.rad),  and  so  on  for  each  parameter,  and  then  the 
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algebraic  solver  is  called  to  solve  the  resulting  equation.  This  can  be  seen 
in  the  put  circle  call  to  the  circle  box.  See  Section  4.1.1  for  a  discussion  of 
the  algebraic  solver  that  the  language  calls  upon. 

6.7     The  JUNO  graphical  editor 

Greg  Nelson's  JUNO  system  is  preceded  by  METAFONT,  the  Xerox  Alto 
Draw  program,  van  Wyk's  IDEAL  language,  and  Borning's  THINGLAB 
system.  The  system  embeds  a  compact  graphical  language  based  on  the 
"guarded  command"  paradigm  of  Dijkstra  [27]. 

6.7.1     Elements  of  the  language 

Some  of  the  elements  of  the  language  are  as  follows: 

Binding.     A  binding  context  is  established  with 

LET  Variables  \  Conetraints  IN  Command  END 

Variables.  Variables  denote  points  in  the  plane.  A  variable  may  be 
introduced  in  a  binding  construct  by  just  giving  its  name,  say  X.  Or  one  may 
give  an  initial  "guess"  u  for  a  point,  which  is  used  for  disambiguation  when 
constraint  satisfaction  results  in  multiple  roots  for  the  system  of  equations 
denoted  by  the  constraints.  This  guess  is  relative  to  a  coordinate  system 
defined  by  a  pair  of  points  {y,z),  where  y  denotes  the  origin,  and  z  the 
endpoint  of  the  unit  vector  (1,0).  This  is  expressed  by 

X==  u  REL  {y,z) 

Geometric  constraints.  Constraints  are  defined  on  line  segments 
denoted  by  pairs  of  points  (i,y)  and  {u,v).  One  can  say  that  two  line 
segments  are  congruent  (same  length)  or  parallel,  and  one  can  say  that  a 
line  is  either  vertical  or  horizontal.  The  syntax  is: 

(i,y)CONG(u,v) 

(i,y)PARA(u,t;) 

HOR(u,t;) 

VER  {u,v) 
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Path  descriptors.  A  path  is  a  sequence  which  may  be  mixture  of 
Hne  segments,  denoted  by  pairs  of  points  {x,y),  and  cubic  Bezier  arcs, 
denoted  by  quadruples  of  points  {x,y,w,z).  Paths  are  used  by  the  graphics 
commands  that  follow. 

Graphics  commands.     The  two  most  basic  drawing  commands 

•  Draw  a  line  along  a  path:  STROKE  Path. 

•  Fill  the  region  circunoscribed  by  a  path:  FILL  Path. 

6.7.2      Juno  code  for  an  equilateral  triangle 

For  example,  the  following  code  will  draw  an  equilateral  triangle  when 
executed  by  the  J  UNO  graphical  editor  [27,  page  239]: 

LET  e  ==  (1.  1)  REL  (a.  b) 

I  (a.  e)  CO^G  (a.  b) 
AND  (b.e)  CONG  (a.  b) 
IN  DRAW  (a.  b).  (b.  e).  (e,  a) 
END 

6.8      The  BERTRAND  production  system  language 

BERTRAND  is  a  production-system  like  language  with  single  assignment  to 
variables  and  types  with  single  supertypes  [21].  Prograjns  consist  of  sets  of 
rules,  together  with  a  special  main  rule  which  is  always  reduced  first. 

6.8.1      Rules  in  BERTRAND 

Rules  are  expressed  as  heads  followed  by  bodies  in  brackets,  with  labels 
and  type  tags. 
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6.8.2     Types  in  BERTRAND 

Types  of  objects  are  either  expressed  by  type  tag  or  by  constructors  created 
by  a  certain  patterned  use  of  rules.  For  example,  to  declare  that  N  is  a  point, 
first  define  points: 

aPoint  {  x:  aNumber  ;  y:  aNumber  ;  true  }  'point 

Then  N  may  be  declared  to  be  a  point  in  one  of  two  ways: 

N:  aPoint 

or 

N'point 

In  the  latter  case,  what  happens  is  that  the  type  point  is  looked  up  in  a 
runtime  table,  and  N'point  is  rewritten  to 

N.x:  aNumber  :  N.y:  aNumber  :  true 

This  is  also  what  happens  in  the  case  N:  aPoint,  except  that  in  this  case,  it 
is  a  consequence  of  the  direct  expansion  of  the  body  of  the  rule  with  head 
aPoint,  without  a  table  lookup  into  the  typve  space. 

BERTRAND  provides  an  assignment  operator,  is.  Expressions  involving 
the  equals  sign  =,  when  they  are  rewritten  such  that  there  is  a  single 
variable  on  one  side,  are  converted  into  the  is  form.  The  next  step  of 
rewriting  stores  the  assignment  into  a  table,  and  replaces  the  assignment 
statement  with  the  value  true.  Every  occurrence  of  the  variable  in  the 
rule  memory  is  then  replaced  by  the  assigned  value,  and  hence  the  assigned 
variable  is,  in  essence,  forgotten,  thus  ensuring  single  assignment  semantics. 
An  example  Is  shown  in  the  following  sequence  of  rewritings,  which  take 
place  in  the  presence  of  a  set  of  rules  for  solving  simple  algebraic  equations: 

x  =  y-|-5;x  =  yx2:y 
xisy  +  5;x  =  y  x2:y 
true  :y  +  5  =  yx2;y 
true  :y-f-5  =  yx2:y 

The  defining  rule  is  that  the  effect  of  an  assignment  of  the  form  x  is  y 
is  the  same  as  adding  the  rule  x  {  y  }  to  the  rule  space. 
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6.8.3      Evaluation  order  in  Bertrand 

The  order  of  evaluation  of  subexpressions  is  undefined  in  BertraND,  in 
order  to  avoid  inducing  a  procedural  semantics.  For  example,  here  is  a 
program  in  BertraND  to  compute  the  factorial  of  S: 

main  {  fact  8  ) 

fact  N  {  prod(ints(N))  } 

ints  1(1} 

ints  A'constant  {  (ints  A-1),  A  } 

prod  A'constant  {  A  } 

prod  (A  .  B)  {  (prod  A)  *  (prod  B)  } 

and  here  is  its  execution  trace,  as  given  by  Leler  (pages  66-67): 

(  prod  (  ints  7-1,7))*  (prod  8) 

(8  *  7)  *  ((  prod  (  ints  6  -  1))  *  (  prod  6)) 

(56  *  6)  *  ((  prod  (  ints  5  -  1))  *  (  prod  5)) 

(336  *  5)  *  ((  prod  (  ints  4  -  1))  *  (  prod  4)) 

(1680  *  5)  *  ((  prod  (  ints  3  -  1))  *  (  prod  3)) 

(6720  *  5)  *  ((  prod  (  ints  2  -  1))  *  (  prod  2)) 

40320 
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7      Conclusion 

The  purpose  of  this  paper  has  been  to  do  a  historical  survey  of  pragmatic 
work  in  the  area  of  graphical  editing  via  the  definition  and  manipulation 
of  constraints  on  geometric  objects.  We  are  seeking  evidence  for  the  ideal, 
namely  a  topologically  correct  editor.  It  seems  to  be  the  case  that  none 
of  the  editors  and  other  graphical  systems  (such  as  typesetting  languages) 
have  had  this  ideal  in  mind.  All  of  these  systems  rely  on  approximate 
nimaerical  computations  which  do  not  support  our  ideal.  However,  we  have 
surveyed  relevant  algorithmic  techniques,  such  as  Grobner  bases,  algebraic 
number  arithmetic,  and  error-preserving  solid  modelling  algorithms,  in  an 
attempt  to  establish  that  the  tools  necessary  for  our  goal  do  exist. 

The  author's  overall  research  objective,  of  which  this  survey  is  a  part, 
is  to  produce  a  working  topologically  correct  graphical  editor.  A  design 
outline  of  this  system  and  the  associated  editing  language  are  in  the  works 
[11].  The  author  hopes  to  complete  a  working  prototype  by  May,  1987, 
and  a  PhD  thesis  on  the  results  by  January,  1988.  Major  stumbling  blocks 
will  be  the  implementation  of  efiBcient  Grobner  bases  and  algebraic  number 
packages,  and  the  implementation  of  "correct  display"  routines  which  do 
not  allow  the  human  viewing  the  display  to  derive  incorrect  information 
about  the  objects  being  projected  onto  the  screen.  This  comes  about  when 
limitations  in  screen  resolution  result  in  incorrect  visual  approximations 
both  on  the  part  of  the  display  algorithm  and  on  the  part  of  the  human 
viewer. 
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